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EXECUTIVE  SUMMARY 


This  final  report  documents  the  results  of  the  research  conducted  at  Advanced  Decision 
Systems  (ADS)  under  a  project  entitled,  “Discrimination  Architecture  Engineering  Support,” 
sponsored  by  Lincoln  Laboratory,  Massachusetts  Institute  of  Technology,  MA .,  for  the  period 
from  December  1,  1987  to  August  31,  1988. 

The  research  relates,  in  general,  to  performance  evaluation  of  systems  for  tracking  and 
classifying  a  large  number  of  dense  (closely-spaced)  objects,  and  in  particular,  to  tracking- 
surveillance  and  object-discrimination  in  a  multi-layer  ballistic  missile  defense  system.  The 
main  purpose  of  the  project  is  to  identify  a  suitable  set  of  performance  measures  and  to 
develop  simple  analytic  models.  It  is  required  that  such  models  be  capable  of  taking  various 
problem  parameters,  such  as 


1.  threat  parameters,  including  threat  trajectories  and  local  densities, 

2.  sensor  parameters,  including  sensor  resolutions,  measurement  accuracies,  detection 
probabilities,  and  false  alarm  rates, 

3.  other  parameters,  including  geometrical  relationships  (viewing  angles,  sensor  coverage, 
etc.),  tracking  time,  number  of  frames,  and  sampling  intervals 


and  converting  them  into  a  suitable  tracking  performance  measure  without  performing  ex¬ 
tensive  Monte-Carlo-type  simulations. 

It  was  decided  that  track  purity  is  the  most  appropriate  tracking  performance  measure 
mainly  because  it  is  a  good  measure  for  representing  the  quality  of  tracking  outputs  especially 
for  object  discrimination  purposes.  This  measure  would  be  important  in  many  tracking 
situations  where  target  classification  is  only  possible  after  correlating  data  from  different 
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time/sensor  samples.  Track  purity  (the  percentage  of  correct  measurements  in  a  given  track) 
can  be  decomposed  into  a  series  of  performance  evaluation  of  sequential  correlation  processes 
(e.g.,  track-to-measurement  associations).  Thus,  a  major  part  of  the  research  effort^  yretc 
devoted  to  establishing  a  simple  analytic  model  for  predicting  the  performance  of  single-scan 
track-to-measurement  correlation.  We  have  derived  and  validated  the  following  scaling  law 
for  predicting  the  probability,  Pc ,  of  a  track  being  associated  with  the  correct  measurement: 

p  _  _  -*4  _  —  v 

rc  =  e  —  c  =  e 

In  this  formula,  the  probability  is  expressed  by  the  object  density  (3  on  the  sensor’s  focal 
plane  and  the  average  prediction  error  standard  deviation  a  (the  error  includes  both  tracking 
(i.e.,  prediction)  error  and  sensor  measurement  error)  or,  equivalently,  the  object  density  /3 
normalized  by  the  error,  or  the  expected  number  v  of  measurements  in  the  one-<r  area  of  the 
error  ellipse. 

In  general,  this  simple  expression  models  the  performance  of  optimal  solutions  to  the 
classical  assignment  problem  in  a  probabilistic  sense.  The  results  of  fairly  extensive  Monte- 
Carlo  simulations  showed  that  this  law  predicts  the  performance  very  well  over  a  reasonably 

large  range  of  assignment  problems.  We  should  note  that  this  evaluation  formula  is  not  meant 

v 

to  be  applied  to  any  specific  sensor  system  and/or  tracking  algorithms.  Rather,  it  expresses 
a  probabilistic  upper  bound  on  the  association  performance  when  only  two  parameters,  i.e., 
object  density  and  track  prediction  errors,  axe  given.  We  have  subsequently  found  out  that 
an  appropriate  expression  when  false  alarms  are  present  is 

pc  —  g— *(0+20jn)9i 

where  /3 fa  is  the  density  of  false  alarms  on  the  sensor’s  focal  plane.  This  model  was  also 
verified  by  Monte-Carlo  simulations. 

With  these  basic  models,  a  method  for  estimating  average  track  purity  given  threat /sensor 
parameters  was  developed.  The  method  depends  on  Cramer-Rao  type  bounds  on  the  track¬ 
ing  accuracies  and  on  dynamical  models  for  predicting  the  object  spatial  densities.  Again 
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the  Cramer-Rao  type  bounds  axe  used  to  obtain  estimated  upper  bounds  on  tracking  per¬ 
formance,  thereby  viewing  potentially  complicated  tracking  processes  (e.g.,  track  initiation, 
continuation,  editing,  elimination,  multiple- hypothesis,  etc.)  as  large  black  boxes  and  mak¬ 
ing  the  resultant  expressions  independent  of  any  particular  tracking  algorithm.  To  verify 
the  track  purity  prediction  method  we  have  developed,  a  small-scale  simulation  of  ballistic 
object  tracking  with  an  orbiting  sensor  was  conducted.  More  extensive  tracking  simulations 
or  examination  of  performance  of  various  tracking  systems  are  planned  for  the  next  phase  of 
this  project. 


Issues  concerning  limited  sensor  resolutions,  i.e.,  closely  spaced  objects,  were  also  inves¬ 
tigated.  Since  our  main  object  is  to  establish  a  method  for  estimating  track  purity,  it  is 
appropriate  to  view  unresolved  measurements  as  bad  measurements  (just  as  incorrectly  as¬ 
sociated  ones).  Thus,  we  have  reduced  the  issues  to  a  problem  of  predicting  the  frequencies 
of  measurements  being  merged  together,  as  functions  of  object  densities  on  focal  planes  and 
sensor  resolutions.  We  have  developed  a  simple  expression  for  the  probability  of  an  object 
being  included  into  unresolved  measurements.  This  simple  expression  was  also  tested  by 
extensive  Monte-Carlo  simulations  and  found  to  be  fairly  accurate. 
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1.  INTRODUCTION 


This  final  technical  report  documents  the  results  of  the  research  which  was  conducted 
at  Advanced  Decision  Systems  (ADS)  under  a  project  entitled,  “Discrimination  Architecture 
Engineering  Support,”  sponsored  by  Lincoln  Laboratory,  Massachusetts  Institute  of  Technol¬ 
ogy,  MA.  The  task  is  to  develop  an  analytical  model  which  relates  scan-to-scan  correlation 
performance  to  threat  and  sensor  parameters.  Since  the  emphasis  is  to  to  develop  analytical 
methods  for  evaluating  sensors’  ability  in  tracking  threat  objects,  the  research  effort  at  ADS 
is  called  “  Tracking  Performance  Evaluation .” 

1.1  Background 

ADS’  task  is  a  part  of  Lincoln  Laboratory’s  effort  on  Discrimination  Architecture  for  the 
Sensors  Directorate  of  SDI  Office.  This  effort  is  concerned  with  modeling  the  interaction  of 
different  sensors  operating  together  in  various  layers  of  a  multi-layer  ballistic  missile  defense 
system.  A  major  element  of  such  a  system  is  a  passive  IR  sensor  operating  in  a  track-while- 
scan  mode  and  covering  a  large  portion  of  the  threat  volume. 

Such  a  system  encounters  various  serious  problems  including: 

•  a  large  number  vof  objects  to  be  tracked  and  discriminated 

•  non-uniform  distribution  of  object  locations  and  velocities 

•  low  observability  -  angle-only  measurements  (no  range  information) 

•  very  high  local  object  density  compared  with  sensors’  resolution 

•  nonlinear  dynamics  (ballistic  trajectories)  and  nonlinear,  time- varying  observation  (or¬ 
biting  sensors) 

•  relatively  long  sensor  revisit  time  compared  with  the  time  constant  of  object  dynamics 

There  is  a  need  to  estimate  the  sensors’  ability  in  tracking  objects  in  this  difficult 
tracking  environment.  For  several  reasons,  it  is  desired  to  perform  such  estimation  without 
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conducting  simulations  using  detailed  threat/sensor  simulators  and  particular  correlation  al¬ 
gorithms.  One  of  the  reasons  is  to  avoid  expensive  Monte  Carlo  runs.  Other  reasons  include 
the  desire  to  estimate  the  upper  bound  of  the  performance  given  basic  threat /sensor  param¬ 
eters  (independent  of  particular  correlation  algorithms),  and  to  produce  sensor  parameter 
specifications  to  meet  given  tracking  performance  requirements. 

Thus  we  need  a  simple  analytic  model  relating  scan-to-scan  correlation  performance  as 
a  function  of  threat  and  sensor  parameters,  which  include: 

•  number  of  objects  to  be  tracked  and  discriminated 

•  local  spatial  density  of  objects 

•  geometric  relationship  between  threats  and  sensors,  such  as  viewing  angles,  detection 
opportunities,  etc. 

•  sensor  resolution,  measurement  accuracy,  detection  probability  and  false  alarm  rate 

•  tracking  time,  number  of  frames,  sensor  revisit  time,  etc. 

Since  we  may  want  to  establish  an  upper  bound  on  tracking  performance,  inclusion  of  op¬ 
erational  parameters  given  to  trackers  (correlators),  such  as  validation  gate  size  (x2-value), 
various  threshold  values,  etc.,  may  or  may  not  be  desirable. 

1.2  Objectives 

ADS’  task  is  to  develop  simple  analytic  models  which  relate  correlation  performance  to 
key  threat/sensor  parameters.  In  particular,  it  is  our  goal  to  develop  scaling  laws  to  identify 
situations  where  similar  performance  can  be  expected  and  to  permit  rapid  comparison  of 
different  scenarios  to  determine  which  will  have  better  scan-to-scan  correlation  performance. 
In  order  to  achieve  this,  it  is  necessary  to  develop  simple  measures  of  scan-to-scan  correlation 
performance.  Analytic  models  and  scaling  laws  thus  obtained  must  then  be  validated  and, 
if  necessary,  calibrated  through  a  series  of  small-scale  simulations. 

Figure  1-1  is  a  conceptual  block  diagram  showing  all  the  system  components,  and  key 
functions  and  parameters.  Our  objective  is  to  establish  simple  analytic  models  which  relate 
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EVALUATION  RESULTS 


Figure  1-1:  Tracking  Performance  Evaluation 


all  the  input  parameters,  going  into  system  components  from  above  in  Figure  1-1,  to  the 
tracking  (correlation)  performance  metrics  coming  out  at  the  bottom  of  the  figure.  One  can 
view  the  whole  system  shown  in  Figure  1-1  as  a  system  which  detects  objects  and  estimated 
their  (six-dimensional)  states.  From  this  point  of  view,  it  is  natural  to  consider  performance 
metrics  such  as 

•  Percentage  of  missing  objects  (or  fraction  of  objects  which  the  system  fails  to  track) 

•  Number  of  false  (extraneous)  tracks 

•  Tracking  accuracy  -  six-dimensional  state  estimation  errors 

which  are  parallel  to  the  standard  set  of  performance  evaluation  metrics  used  for  sensors, 
i.e.,  detection  probability,  false  alarm  rate  and  measurement  accuracy.  Unlike  sensor  signal 
processors  which  usually  operate  on  a  snap-shot  (single-scan)  basis,  trackers  or  correlators 
involve  very  complicated  data  processing,  and  as  a  result,  it  is  sometimes  not  straightforward 
as  to  how  to  define  the  above  three  metrics. 

In  our  particular  environment,  i.e.,  in  ballistic  missile  defense  systems,  however,  it  is 
more  appropriate  to  measure  the  third  criterion,  i.e.,  tracking  accuracy,  as  the  closeness 
of  the  data-to-data  correlation  (partitioning  of  multi-scan  detections)  results  to  the  real 
source  of  each  measurement.  This  is  so  because,  in  the  present  stage  of  technology,  effective 
object  discrimination  requires  time  series  of  measurements  distinctively  originating  from 
each  object  to  be  discriminated.  In  other  words,  discrimination  algorithms  must  rely  on  the 
correlation  results  of  object  trackers.  In  such  a  situation,  mis-association  generated  by  a 
tracker  (correlator)  is  naturally  a  serious  problem  for  a  discrimination  system,  and  it  is  very 
important  to  predict  the  magnitude  of  such  mis-association.  Thus  our  main  objective  is  to 
predict  tracking  performance  in  terms  of  how  pure  each  track  is,  in  other  words,  how  many 
mis-associations  each  track  contains. 

1.3  Technical  Approach 

As  mentioned  before,  the  primary  objective  of  this  research  is  to  quantify  the  (poten¬ 
tial)  ability  of  tracking  and  discriminating  objects  given  a  set  of  threat/sensor  parameters. 
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Therefore,  it  is  desirable  to  conduct  analysis  independent  of  any  specific  tracking  (correlation) 
algorithm.  One  way  to  accomplish  this  objective  is  to  analytically  evaluate  the  performance 
of  the  optimal  algorithm.  A  class  of  general  multi-target  tracking  problems  has  been  formu¬ 
lated  and  optimal  solutions  to  them  have  been  developed.  (See,  e.g.,  [l]).  Unfortunately, 
performance  evaluation  based  on  this  general  theory  of  multi-target  tracking  is  extremely 
difficult  even  in  Gauss-Poisson  cases.  This  is  so  because  such  analysis  cannot  avoid  com¬ 
plicated  probabilistic  calculation  involving  non-gaussian  (typically  x2)  random  variables.  In 
fact,  almost  all  the  tracking  algorithms  (including  the  tracking  algorithm  specifically  de¬ 
signed  to  solve  the  tracking  problem  in  ballistic  missile  defense  systems  [2])  involve  x2-type 
random  variables  in  one  way  or  another,  which  makes  analytic  evaluation  very  difficult. 

Nonetheless,  there  have  been  a  small  number  of  papers  and  reports  on  this  subject 
mainly  considering  specific  algorithms,  e.g.,  [3],  [4]  and  [5].  Primary  concerns  were  on  statis¬ 
tics  of  missing  objects  and  false  tracks ,  and  analytic  results  were  obtained  from  Markov  chain 
models  which  approximately  describe  detection  history  for  each  track.  Therefore,  emphasis 
was  placed  on  track  initiation  performance  rather  than  track  continuation ,  and  the  possibil¬ 
ity  of  mis-association  was  not  considered.  In  the  context  of  ballistic  missile  defense  systems, 
or  at  least  of  in  some  phases  of  such  systems,  failure  of  track  initiation  or  creation  of  (perma¬ 
nent)  false  tracks  are  not  of  central  importance  due  to  relatively  high  detection  probability 
and  substantially  regular  dynamics.  Nonetheless,  the  track  initiation  itself  is  quite  important 
due  to  the  high  object  density,  the  limited  resolution  and  the  low  observability  (angle-only 
measurements).  On  the  other  hand,  the  track  initiation  problem  may  be  made  somewhat 
easier  through  system  architecture  concepts  such  as  birth-to-death  tracking.  On  the  other 
hand,  track  quality  in  terms  of  track  purity  is  of  great  importance  for  object  discrimination 
and  targeting  purposes,  and  is  thus  the  performance  measure  used  in  this  study. 

Since  track  purity  may  be  expressed  as  a  summation  of  the  correctness  of  track-to- 
measurement  correlation  at  each  scan,  the  most  critical  analysis  is  that  on  the  single-scan 
correlation  performance.  Our  approach  is  first  to  establish  a  very  simple  analytic  model 
which  relates  the  single-scan  correlation  performance  to  a  few  parameters  characterizing  the 
scan  in  question,  and  then  to  relate  that  single-scan  model  to  prediction  of  average  track 
purity.  Figure  1-2  summarizes  our  technical  approach  for  this  research.  In  order  to  develop 


a  simple  analytic  model,  we  must  base  the  single-scan  analysis  on  a  relatively  simple  math¬ 
ematical  model  with  several  simplifying  assumptions.  However,  crucial  parameters  must  be 
adequately  included.  In  our  approach,  the  single-scan  correlation  performance  is  related  to 
the  measurement  prediction  accuracy  by  tracks  and  object  density  on  the  sensor’s  focal  plane 
in  each  scan.  These  two  key  parameters  can  then  be  related  to  other  parameters  defined  by 
threat/sensor  scenarios.  In  arriving  at  analytic  expressions  for  correlation  performance,  it 
is  inevitable  to  simplify  models  by  ignoring  complicated  factors  and  approximating  math¬ 
ematical  formulae.  In  general,  in  our  approach,  initially  analysis  is  conducted  with  many 
assumptions  (sometimes  restrictive),  and  then,  relaxation  of  assumptions  is  attempted  later 
one  by  one.  Each  step  of  analytic  results  will  be  verified  as  much  as  possible  through  Monte 
Carlo  simulations. 

The  next  section,  Section  2,  describes  the  single-scan  correlation  performance  analy¬ 
sis  which  forms  the  basis  of  the  majority  of  the  results  described  in  this  report.  The  next 
section,  Section  3,  shows  how  key  parameters  used  in  Section  2  can  be  derived  from  a 
given  threat /sensor  scenario;  this  is  followed  by  the  expressions  for  track  purity  in  Section 
4.  Throughout  these  sections,  various  kinds  of  Monte  Carlo  simulations  are  used  to  verify 
the  analytic  expressions.  The  last  section,  Section  5,  summarizes  the  results  and  discusses 
directions  for  future  research.  Detailed  mathematical  derivations  are  presented  in  the  Ap¬ 
pendices. 


7 


2.  SINGLE-SCAN  ANALYSIS 


This  section  describes  the  single-scan  correlation  analysis  which  serves  as  the  basis  of  the 
studies  presented  in  this  report.  In  a  sense,  this  section  describes  the  problem  of  estimating 
track  purity  in  its  simplest  form.  Initially,  we  will  ignore  some  factors,  i.e.,  missed  detections, 
resolution  limitation  (closely  spaced  objects),  false  alarms,  etc.  The  effects  of  those  factors 
will  be  considered  later  in  this  section  or  in  the  subsequent  sections.  The  mathematical  model 
used  in  this  section  is  almost  identical  to  that  in  [6],  which  analyzes  the  single-scan-based 
sensor-to-sensor  correlation  problem.  The  technical  report  [6],  written  fourteen  years  ago  is 
concerned  with  the  precise  calculation  of  the  probability  of  the  most  likely  correlation  being 
the  correct  one  and  draws  its  main  results  through  the  Fourier  (characteristic  function) 
technique.  In  this  section,  by  introducing  further  simplifying  assumptions,  a  very  simple 
expression  of  the  probability  will  be  obtained  as  a  basic  scaling  law. 

2.1  Mathematical  Model  for  Single-Scan  Correlation  Problem 

Suppose  that  we  have  been  tracking  N  objects  in  a  number  of  scans  from  a  sensor 
(or  sensors)  from  which  N  tracks  have  been  formed,  and  that  there  is  a  perfect  one-to-one 
matching  between  the  set  of  “real”  objects  and  that  of  tracks.  We  assume,  for  the  moment, 
there  are  neither  false  alarms  nor  missed  detections  in  the  previous  scams  as  well  as  the  scan  to 
be  considered  here.  Thus  we  consider  a  sensor  scan  consisting  of  N  measurements,  yi, 
with  each  y,  originating  from  the  i-th  object.1  Assuming  an  angle-only  measurement  passive 
sensor,  each  y,  is  modeled  as  a  two-dimensional  random  vector, 

y>  =  z{  +  n{  (2.1) 

for  each  i  £  {1,...,1V},  where  z;  is  the  “true”  projection  of  the  position  of  the  t-th  object 
onto  the  sensor’s  focal  plane  and  n{  is  the  additive  zero-mean  gaussian  noise  with  a  2  x  2 
variance  matrix  ZZ,. 

lIn  an  actual  scan,  it  is  natural  to  model  the  order  of  measurements  as  a  realisation  of  a  totally  random 
(i.e.,  equally  probable)  permutation  of  a  given  object  labeling.  Here  we  are  using  the  identity  index  map 
from  the  set  of  the  objects  to  the  set  of  measurements  for  notational  convenience.  For  a  correlator,  the 
identity  map  is  just  one  of  all  the  equally  probable  permutations.  See  Appendix  A  for  a  related  discussion. 
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On  the  other  hand,  for  each  i  €  {1,...,JV})  the  i-th  irack  (tracking  the  i-th  object)  is 
associated  with  measurement  prediction  y,-  which  is  modeled  as 

Vi  =  Zi  +  mi  (2.2) 

where  m,-  is  a  zero-mean  gaussian  random  vector  representing  the  measurement  prediction 
error  of  the  i-th  track,  with  Qi  being  its  variance  matrix.2  We  assume  that  the  random 
vectors,  yi ni ,...,njv,  are  independent  from  each  other.3 

In  this  simplified  situation,  it  is  natural  to  evaluate  the  performance  of  the  optimal 
assignment  algorithm.  Each  assignment  is  a  permutation  7r  on  the  index  set  {1, ...,  jY}, 
i.e.,  under  a  permutation  x,  the  i-th  track  is  assigned  to  the  x(i)-th  measurement.  Under 
an  appropriate  set  of  assumptions,  it  can  be  shown4  that  an  optimal  track-to-measurement 
assignment  is  represented  by  a  permutation  which  maximizes  the  assignment  cost  function 

s  ' 

A*)  =  £x(*»  *(*))*  (2-3) 

»=i 

over  the  set  of  all  the  permutations  x  on  the  set  {1, N},  where,  for  each  (i,  j), 

=  lift  -  (2-4) 

with5 

=  Qi  +  Ri  (2-5) 

which  is  assumed  to  be  positive  definite.  Then  the  probability  of  the  optimal  assignment 
being  the  “true”  one  is  the  probability  of  the  event  in  which  J(t)  <  J(ir)  for  all  the  permu¬ 
tations  7r  where  i  is  the  identity  map6  representing  the  “true”  assignment. 

The  problem  of  determining  an  optimal  7r  to  minimize  (2.3)  is  the  well-known  classical 
assignment  problem.  It  is  also  a  well  known  result  that,  under  a  correct  assignment  (x  =  i), 

2  In  an  actual  situation  where  ballistic  objects  are  being  tracked,  each  measurement  prediction  in  and 
its  associating  error  variance  matrix  Qi  are  obtained  by  the  nonlinear  projection  from  the  object  state 
distribution  for  each  track  on  to  the  sensor’s  focal  plane. 

3See  Appendix  A  for  discussions  concerning  this  assumption. 

4See  Appendix  A. 

5 By  jjzjj.4,  we  mean  a  semi-norm  of  z  in  a  Euclidean  space  defined  by  a  nonnegative  definite  symmetric 
matrix  A  as  Ijzljx  d=  \/zT Az.  Thus,  when  I  is  the  identity  matrix,  ||z||  =  |jz||/  is  the  Euclidean  norm  of  z. 
zT  is  the  transpose  of  a  vector  or  matrix  z. 

•Namely,  *.(«)  =  t  for  all  t  €  {1, ..,  N} 
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the  assignment  cost  J(l)  is  a  x2  random  variable  with  degree  of  freedom  equal  to  N  times 
the  dimension  (two)  of  the  measurement  space.  It  can  also  be  shown  that,  when  ir  is  not  a 
true  association,  J(ir)  is  a  non-centric  x2  random  variable  (See  [6].).  Although  these  results 
have  been  known  among  researchers  for  at  least  two  decades,  there  is  not  much  else  in  the 
form  of  analytical  facts  that  we  can  point  out. 

2.2  Estimation  of  Probability  of  Correct  Association 

In  [6],  it  was  argued  that  the  events  in  which  three-way  or  more  complex  assignment 
errors  happen  in  an  optimal  assignment  are  very  rare  so  that  we  should  only  consider  two- 
way  mis-associations.  This  argument  may  be  loosely  supported  by  the  following  observation: 
Every  permutation  can  be  expressed  as  a  composition  of  transpositions  (two-way  switches). 
Therefore,  unless  the  probability  of  having  each  transposition  is  significantly  high,  two  or 
more  simultaneous  transpositions  involving  the  same  object  should  hot  happen  very  often. 
It  should  be  very  rare  particularly  when  the  object  density  is  not  significantly  high.  For  the 
moment,  let  us  assume  that  transpositions  are  the  only  significant  events. 

Single  Transposition:  To  simplify  the  discussion,  let  us  consider  each  transposition 
independently,  and  thus,  assume  only  two  objects,  i.e.,  N  —  2.  The  event  to  consider  is 

Eu  =  |j({(l,l),(2,2)})  <  J({(1,2),(2,1)»  j  (2.6) 

of  two  tracks  being  correctly  associated  with  measurements.  As  seen  in  [6]  and  in  Appendix 
B,  the  calculation  of  this  probability  is  not  easy  and  evaluation  may  be  possible  only  through 
numerical  methods,  e.g.,  Fourier  transformation  (i.e.,  by  evaluating  the  characteristic  func¬ 
tion  of  a  non-centric  x2  random  variable),  which  will  not  lead  to  any  analytic  expression 
useful  for  our  purpose.  However,  when  we  assume  Qj  =  Qi  and  R\  =  /?2,  and  when  condi¬ 
tioned  by  the  predictions  by  tracks,  we  have  a  simple  probabilistic  expression, 

=  erf  ^ 

where7 

SijSS  (2.8) 

7erf()  is  the  error  function  defined  as  erf(*) 


Prob(E12|t/i,j/2) 


v/21 


/i  ~  Sills-1 


(2.7) 


10 


Figure  2-1:  Single-Scan  Analysis:  Two  Objects 

for  all  (t,j)  €  {1,2}  x  {1,2}.  This  assumption,  (2.8),  may  be  justified  in  situations  where 
object  states  in  a  group  are  very  similar  and  so  are  the  SNR’s  of  each  detection.  It  may  not  be 
a  very  plausible  assumption,  however,  in  special  situations,  e.g.,  when  two  object  groups  are 
crossing  or  merging  on  a  sensor’s  focal  plane.  Even  in  such  situations,  with  enough  number 
of  frames  being  accumulated,  this  assumption  may  become  eventually  a  close  approximation 
of  the  reality. 

Assume  further  that  yi  —  y2  is  distributed  uniformly  on  a  disk  with  radius  r,  or  equiv¬ 
alently  that  yi  is  at  the  origin  0  and  y%  is  distributed  uniformly  on  the  disk  with  radius  r 
(See  Figure  2-1.).  Define  an  average  prediction  error  standard  deviation  cr  as  the  geometric 
average  of  the  semi-major  and  the  semi-minor  of  the  ellipse  defined  by  jjy||s-i  =  1,  i.e.8, 

er  =  [det(S)]1/4  (2.9) 

Then,  as  shown  in  Appendix  B,  when  the  radius  r  is  large  enough  compared  with  the  square 
root  of  maximum  eigenvalue,  i.e.,  the  semi-major  of  the  ellipse  defined  by  5,  the  unconditional 
*By  det(/i),  we  mean  the  determinant  of  a  square  matrix  A. 
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Figure  2-2:  Single-Scan  Analysis:  N  Objects 


probability  of  the  event  Eu  can  be  expressed  simply  as 

Prob.(f^2)  «  1  -  ^  *  o  (2.10) 

Multiple  Transpositions:  Now  let  us  go  back  to  our  iV-object  model.  Without  loss 
of  generality,  we  can  place  the  track  we  are  concerned  with,  say  the  first  track,  at  the  origin, 
i.e.,  yi  =  0.  Assume  that  all  the  other  tracks  are  uniformly  distributed  on  the  disk  with 
radius  r  (See  Figure  2-2).  This  assumption  should  not  be  too  restrictive  as  far  as  the  radius 
r  is  large  enough  compared  with  the  average  prediction  error  standard  deviation  a  and  the 
track  we  are  considering,  i.e.,  jh,  is  not  almost  on  the  edge  of  the  track  group’s  extent. 
Furthermore,  let  us  assume  (2.8)  for  all  pairs  (t,j)  of  indices  for  tracks  and  measurements. 
Let  Ee  be  the  event  in  which  there  is  no  transposition  between  the  first  track  t/x  and  any  of 
the  remaining  tracks.  Then,  by  assuming  that  any  mis-assodation  is  caused  by  one  or  more 
such  transpositions  involving  the  first  track,  the  probability  of  the  first  track  being  correctly 
associated  with  the  measurement  can  be  simply  expressed  (following  from  (2.10))  sis 

Prob.(Ee)  w  (Prob.^ij))*'1  *  ^1  -  «  1  -  *(3c2  (2.11) 
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where  0  is  the  object  density  on  the  focal  plane,  i.e., 


N  -  1  _  N_ 

7Tr2  TTT2 


(2.12) 


assuming  N  1. 

It  follows  from  (2.11)  that,  when  the  number  N  of  objects  is  a  random  integer  with  a 
Poisson  distribution,  the  probability  of  the  track  being  associated  with  the  correct  measure¬ 
ment  becomes 

Prob.(Ee)  Rs  exp  (—ir0 a7^j  (2.13) 

where  the  object  density  0  is  now  re-defined  as 


(2.14) 


with  v  being  the  expectation  of  the  number  N  of  objects,  which  is  now  re-defined  as  a  Poisson 

\ 

random  variable. 

The  expression  (2.13)  is  the  most  important  result  in  the  studies  described  in  this  report. 
It  is  a  simple  scaling  law  which  expresses  the  probability  of  correct  association  as  the  function 
of  1)  object  density  0  on  the  focal  plane  and  2)  the  average  measurement  prediction  error 
standard  deviation  a.  Since  (2.13)  is  a  function  of  0a7,  we  can  have  a  simpler  expression 

Prob.(Ee) «  exp(-*0)  (2.15) 

where  0  is  the  normalized  object  density  defined  by 

0  =  0a7  =  0yJdet(S)  (2.16) 

It  should  be  noted  that  (2.15)  claims  that  the  correlation  performance  is  determined  solely 
by  the  normalized  object  density  (defined  by  (2.16))  on  the  focal  plane. 

Monte  Carlo  Simulation:  To  verify  the  scaling  law  (2.13)  or  (2.15),  a  series  of  Monte 
Carlo  simulations  were  performed.  Figure  2-3  shows  the  results.  For  each  run,  a  hundred 
tracks  were  simulated  by  generating  (y«),-^  uniformly  on  a  disk  with  an  appropriate  size, 
and  then  one  hundred  measurements,  (y,)}^,  were  generated  according  to  eqns.  (2.1)  and 
(2.2).  Then  the  optimal  assignment  was  calculated  by  Munkres  algorithm  ([7]),  and  it  is 
subsequently  compared  with  the  true  track- to- measurement  correspondence.  The  probability 
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Figure  2-3:  Single-Scan  Correlation  Performance  (1) 

of  correct  association  is  defined  as  the  percentage  of  tracks  which  were  correctly  assigned 
with  measurements  and  is  averaged  over  50  samples  for  each  point  in  Figure  2-3.  The 
measurement  prediction  error  variance  matrix  5  was  chosen  to  be  a  diagonal  matrix  with 
the  same  diagonal  elements.  The  cases  with  non-diagonal  matrices,  or  in  other  words,  the 
non-zero  eccentricity  of  the  ellipse  defined  by  [|y||s-«  =  1,  will  be  discussed  later  in  Section 
2.3. 

When  we  derived  the  scaling  law,  (2.13),  we  used  a  mathematical  model  which  places 
the  track  to  be  considered  at  the  origin  and  assumed  other  tracks  are  distributed  uniformly 
around  it.  When  a  track  is  situated  at  the  peripheral  of  a  track  group,  the  track  density 
which  affects  the  probability  of  correct  association  effectively  decreases.  This  kind  of  “edge” 
effect  was  not  considered  in  (2.13).  Since  the  purpose  of  the  Monte  Carlo  simulation  is  to 
verify  the  appropriateness  of  the  simple  mathematical  expression  (2.13),  we  considered  only 
the  tracks  which  are  within  a  disk  with  half  the  radius  of  that  of  the  disk  over  which  all  the 
tracks  are  uniformly  distributed,  thereby  eliminating  the  edge  effects.  We  will  discuss  issues 
concerning  uneven  object  density  later  in  this  section. 
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Figure  2-3  shows  the  results  with  three  different  object  densities,  i.e.,  /3  being  .1,  .25 
and  .5  objects  per  unit  area  (distinguished  by  three  different  symbols  in  the  figure),  as  plots 
of  the  correct  association  probability  versus  the  average  measurement  prediction  variance  <r2 
(measured  by  unit  area),  together  with  the  scaling  law  (2.13).  Figure  2-4  plots  the  same  data 
using  the  normalized  object  density  defined  by  (2.16)  and  compares  them  with  the  scaling 
law  (2.15).  Figure  2-4  uses  the  same  symbols  used  in  Figure  2-3  for  different  actual  (not 
normalized)  object  densities.  In  Figures  2-3  and  2-4,  the  empirical  results  agree  with  our  an¬ 
alytic  results,  eqns.  (2.13)  and  (2.15),  very  well  except  for  regions  where  we  have  high  object 
densities  or  large  measurement  prediction  errors.  In  such  regions,  however,  both  actual  and 
predicted  correlation  performance  are  very  bad,  i.e.,  with  probability  of  correct  association 
being  50%  or  less,  and  hence,  it  is  not  reasonable  even  to  try  track-to-measurement  correla¬ 
tion  at  all.  Nonetheless,  it  is  interesting  to  see  that  the  actual  optimal  assignment  performs 
much  better  than  the  theoretically  expected  performance.  This  deviation  seems  to  originate 
from  the  fact  that,  in  our  analysis,  the  mis-association  mechanism  is  modeled  simply  by 
composition  of  independent  transpositions,  and  complicated  multiple-way  mis- associations 
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Figure  2-5:  Single-Scan  Correlation  Performance  (3) 


and  dependence  of  association  errors  (ignored  in  the  analysis)  cannot  be  ignored  in  high 
object  density  or  large  measurement  prediction  error  cases.  Another  explanation  for  this 
deviation  is  simply  the  increase  of  “edge”  effects  which  we  have  discussed  before.  (See  also 
Appendix  B  for  related  discussions.)  Figure  2-5  shows  the  same  plots  as  those  in  Figure  2-4 
but  the  abscissa  is  taken  as  the  expected  number  of  measurements  in  the  validation  gate  of 
each  track.  The  expected  number  vq  of  measurements  per  track  validation  gate  is  calculated 


as 

vg  =  1  +  */?Xc\/det(S)  (2-17) 


where  xc  is  the  X2'^0®  to  determine  the  validation  gate  size.  In  this  plot,  the  value,  xg  =  3, 
is  used. 

It  may  be  interesting  to  see  that,  in  actual  simulations,  how  often  multiple  switches, 
i.e.,  mis-association  not  caused  by  simple  transpositions,  happen.  Table  2-1  shows  the  av¬ 
erage  number  of  occurrence  of  various  kinds  of  mis-associations  with  different  normalized 
object  densities,  (using  unit  l/x).  As  seen  in  the  table,  the  occurrence  of  complicated 
object  switches  increases  as  the  object  density  increases.  However,  within  the  range  where 
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Table  2-1:  Mis-Association  Statistics 


KBB1 

Two-Way  Switches 

Three-Way  Switches 

Four  Ways  or  More 

.02 

.7 

0 

0 

.05 

1.9 

.1 

0 

.1 

3.8 

.4 

.1 

.2 

6 

1.2 

.1 

.5 

9 

2.8 

1.4 

1 

10.1 

4.7 

4.6 

2 

7.2 

3.6 

6.8 

the  probability  of  correct  association  is  meaningful,  say  about  50  %,  the  two-way  switches 
dominates  more  complicated  switches.  This  is  consistent  with  the  fact  that  the  scaling  law, 
(2.13)  or  (2.15),  considering  only  two-way  switches,  predicts  the  assignment  performance  rel¬ 
atively  well  (as  seen  in  Figures  2-3  through  2-5).  The  effects  of  complicated  switches  which 
are  not  modeled  in  our  analysis  are  not  clear.  The  fact  that  the  simple  expression,  (2.13) 
or  (2.15),  fits  well  with  simulation  results  indicates  overall  cancellation  of  over-estimation  of 
the  performance  with  under-estimation. 

2.3  Effects  of  Eccentricity  of  Error  Ellipse 

In  the  single-scan  analysis  in  this  section,  we  have  assumed  the  prediction  error  variance 
matrices,  5^,  for  all  the  track-measurement  pairs  are  identically  equal  to  a  matrix  5.  The 
matrix  5  (which  is  assumed  to  be  positive  definite)  defines  an  ellipse®,  {y  €  ^.J|||l/|ls-j  =  1}- 
In  the  simulations,  the  results  of  which  are  shown  in  the  previous  section,  2.2,  we  have  as¬ 
sumed  zero  eccentricity.  In  this  subsection,  we  will  examine  the  effects  of  eccentricity  of  the 
error  ellipses  and  measure  such  effects  by  means  of  Monte  Carlo  simulations  to  confirm  the 
validity  of  our  analytic  expression,  i.e.,  (2.13).  According  to  (2.13)  together  with  the  defi¬ 
nition  (2.9)  of  the  average  error  standard  deviation  a,  the  probability  of  correct  association 
can  be  expressed  as 

Prob.(£e)  ft:  exp(—  (5A) 


*72  =  (—  oo,  oo)  =  set  of  real  numbers. 
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[Probability  of  Correct  Association 


Figure  2-6:  Effects  of  Eccentricity  of  Prediction  Error  Ellipses 


where  A  is  the  area  of  the  error  ellipse  defined  by  {y  6  =  1}-  In  other  words,  the 

probability  of  correct  association  can  be  expressed  as 

Prob.(£c)  ss  exp  (Expected  Number  of  Objects  in  Error  Ellipse) 

Let  hi  and  H2  be  the  two  (strictly  positive)  eigen-values  of  the  matrix  5,  and  let  a  = 
m&x{y/Hi,  y/m}  and  b  =  min{^/pT,  yffii}  be  the  semi-major  and  the  semi-minor  of  the 
ellipse.  Then  the  average  prediction  error  standard  deviation  c  defined  by  (2.9)  can  be 
rewritten  as 

a  =  aVl  —  e2  (2.18) 

and  hence,  (2.13)  becomes 

Prob.(jFe)  ~  exp(— 7r/5aJ(l  —  eJ))  (2.19) 

where  e  is  the  eccentricity  of  the  ellipse. 

Figure  2-6  shows  the  probability  of  correct  association  where  the  semi-major  a  of  the 
error  ellipses  and  the  object  density  0  are  fixed  such  that 

0  =  0a2  =  0  a2  =  .05 
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when  the  eccentricity  is  zero  (e  =  0)  and  the  eccentricity  is  varied  from  0  to  1.  This  figure 
also  shows  the  results  of  Monte  Carlo  simulations.  For  each  Monte  Carlo  simulation  point 
shown  in  the  figure,  200  samples  are  taken  for  the  simulation  using  100  objects.  As  seen 
in  the  figure,  the  simulation  results  closely  follow  the  theoretical  curve  shown  by  the  solid 
curve. 

2.4  Uneven  Local  Object  Density 

In  Section  2.2,  we  have  derived  simple  scaling  laws  to  determine  the  probability  of  correct 
association  as  a  function  of  object  density  and  track  measurement  prediction  errors.  To  do 
so,  we  have  assumed  objects  are  distributed  uniformly  on  a  given  region  of  the  sensor’s  focal 
plane.  This  uniformity  assumption  may  not  validly  reflect  reality.  For  example,  if  objects  are 
distributed  uniformly  within  a  ball  in  a  three-dimensional  space,  apparently  their  projections 
onto  a  plane  are  not  distributed  uniformly.  In  general,  the  two-dimensional  projection  of 
a  group  of  three-dimensional  objects  is  such  that  object  density  is  in  general  a  decreasing 
function  of  the  distance  from  the  center  of  the  image.  As  will  be  shown  in  Section  4,  even  in 
such  a  case,  if  the  object  density  is  appropriately  calculated,  the  scaling  law,  (2.13)  or  (2.15), 
still  provides  us  with  relatively  accurate  performance  prediction.  However,  one  may  want  to 
question  how  uneven  local  object  density  affects  the  correlation  performance,  in  particular, 
one  may  want  to  quantify  the  difference  in  correlation  performance  between  inner  objects 
and  outer  ones,  within  the  same  object  group. 

A  naive  way  to  derive  an  expression  which  may  handle  the  dependency  on  the  varying 
object  density  would  be 

Prob.(iTe)  ~  exp  7r/3 {p)&7^j  (2.20) 

or 

Prob. (2?c)  «  exp(— 7 *0{p))  (2-21) 

which  are  obtained  directly  from  (2.13)  and  (2.15)  (resp.)  by  simply  replacing  the  constant 
object  density  0  (normalized  density  0)  by  non-constant  density,  e.g.,  the  density,  0(p)  or 
0(p),  depending  on  the  distance  p  from  the  center  of  the  track  group.  However,  there  is  no 
immediate  justification  to  do  so.  This  is  so  because,  although  we  can  model  varying  object 
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Figure  2-7:  Uneven  Probability  of  Correct  Association 

density  by  a  non-constant  function  /3  or  /3,  such  notion  is  only  valid  for  a  small  region  where 
we  cannot  ignore  the  extent  of  the  prediction  error  ellipse.  Nonetheless,  when  the  object 
density  is  a  monotonic  (decreasing)  function  of  the  distance  from  the  center,  the  effects  of 
the  slope  of  the  object  density  may  be  averaged  out  so  that  expression  (2.20)  or  (2.21)  is  in 
fact  a  good  approximation  for  the  non-uniform  local  object  density. 

In  order  to  illustrate  this  point,  the  probability  of  correct  association  was  calculated 
assuming  a  gaussian  density,  i.e., 

(2-22) 

with  S  <r,  where  p  is  the  distance  from  the  center  of  a  track  group  and  v  is  the  expected 
total  number  of  tracks  in  the  group,  rather  than  assuming  a  uniform  density.  Appendix 
C  shows  (2.20)  is  in  fact  a  good  approximate  expression  in  such  a  case.  This  analysis  is 
also  supported  by  a  Monte  Carlo  simulation,  the  results  of  which  are  show  in  Figure  2-7. 
The  curve  shown  in  this  figure  is  obtained  from  (2.20)  and  (2.22)  with  v  —  20,  S  =  2  and 
a  —  .2.  This  theoretical  curve  is  supported  by  a  Monte  Carlo  simulation  whose  results  are 


20 


shown  in  the  figure  by  step  functions.  The  simulation  produces  300  samples  of  a  20-object 
scenario  according  to  (2.22)  with  parameters  consistent  with  the  curve.  Then  all  the  correct 
associations  and  mis-associations  are  sorted  according  to  the  distance  of  the  track  from 
the  center  of  the  group,  and  the  average  probability  of  correct  association  is  calculated  for 
each  interval  of  the  distance.  The  figure  shows  that  the  simulation  results  coincide  with  the 
theoretical  curve  surprisingly  well,  suggesting  that  this  kind  of  “gaussian”  model  may  be  used 
to  investigate  the  unevenness  within  an  object  group  in  terms  of  correlation  performance. 

2.5  Effects  of  False  Alarms 

In  section  2.2,  we  have  derived  a  very  simple  but  relatively  accurate  expression  for  the 
probability  of  correct  association  by  first  1)  calculating  the  probability  of  transposition  (two- 
way  switch)  of  objects  (assuming  only  two  objects)  and  then  2)  identifying  the  event  of  correct 
association  as  the  one  in  which  no  transposition  takes  place  with  any  of  the  surrounding 
objects.  With  this  approach,  it  is  rather  straightforward  to  take  false  alarms  into  account. 
Recalling  that  Ec  is  the  event  of  a  track  being  associated  with  a  correct  measurement,  let 
Ed  be  the  event  that  there  is  no  transposition  between  the  track  in  question  and  any  other 
surrounding  track  and  Ecf  be  the  event  that  there  is  no  transposition  between  the  track  and 
any  of  the  false  alarms  which  surrounds  it.  Then  we  have 

Prob.(Ee)  =  Prob^EctJProb^Ee/)  (2.23) 

which  may  be  supported  by  the  assumption  that  tracks  (i.e.,  objects)  are  independent  of 
false  alarms. 

For  the  sake  of  simplicity,  let  us  assume  that  the  common  prediction  error  variance 
matrix  5  is  an  identical  matrix  with  the  same  diagonal  element  a2.  Then  as  shown  in 
Appendix  D,  we  have 

Prob. (.£<.*)  «  exp(— 7r /3a2)  =  exp(— 7r/j)  (2.24) 

and 

Prob.(£'e/)  %  exp(—2irfiFA<r7)  (2.25) 

where  (3fa  is  the  density  of  false  alarms.  It  is  interesting  to  observe  that  the  two  expressions 
(2.24)  and  (2.25)  are  almost  identical  with  the  difference  of  the  factor  2,  which  originates 
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Figure  2-8:  Effects  of  False  Alarms 
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from  the  fact  that  measurements  from  other  tracks  contain  their  prediction  errors  while  false 
alarms  are  simply  assumed  to  distributed  uniformly  around  the  track  in  question. 

It  follows  from  (2.23),  (2.24)  and  (2.25)  that 

Prob.(£e)  r;  exp(— t(/3  +  2f3FA)v7)  =  exp(-ir/?(l-f  27Fit)a2)  =  exp(-ir0(l  +  27Fx))  (2.26) 

where  7 fa  is  the  ratio  of  the  false  alarm  density  over  the  object  density  on  the  focal  plane, 
i.e.,  PfaIP-  Figure  2-8  shows  the  effects  of  the  false  alarms  on  the  probability  of  correct 
association.  In  this  figure,  the  theoretical  curve  defined  by  (2.26)  is  shown  as  well  as  Monte 
Carlo  simulation  results  (shown  by  dots).  The  normalized  object  density  /?  of  .1  is  used  for 
this  figure  and  the  ratio  of  the  false  alarm  density  over  the  object  density  is  chosen  as  the 
independent  parameter.  For  the  simulation,  100  targets  are  generated  for  each  run  and  each 
point  is  obtained  form  50  samples.  As  seen  in  this  figure,  simulation  results  support  the 
simple  expression  (2.26).  As  observed  in  Figure  2-3  and  Figure  2-4,  the  theoretical  expres¬ 
sions  give  us  slightly  pessimistic  predictions  when  the  correlation  performance  is  degraded 
substantially  due  to  dense  false  alarms. 
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2.0  Summary  of  Section  2 


In  this  section,  we  have  presented  the  major  results  of  the  studies  described  in  this 
report,  i.e.,  the  scaling  law  expressed  by  (2.13)  and  (2.15).  This  analytical  expression  is 
very  simple  and  yet  we  have  shown,  by  means  of  extensive  Monte  Carlo  simulations,  that 
it  may  provide  very  accurate  prediction  for  single-scan  correlation  performance.  A  crucial 
assumption  in  deriving  this  expression  is  that  the  prediction  error  variance  matrix  for 
each  track-measurement  pair  is  (at  least  almost)  identical,  which  may  be  valid  when  we  track 
a  group  of  objects  whose  state  vectors  are  close  together.  Issues  concerning  the  shape  of  the 
error  ellipses,  the  uneven  local  object  density  and  the  false  alarms,  have  been  also  discussed 
in  this  section.  Issues  concerning  the  unresolved  measurements  (i.e.,  closely  spaced  objects) 
will  be  included  in  Section  4. 

\ 
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3.  TRACK  ACCURACY  AND  OBJECT  DENSITY 

In  Section  2,  we  have  developed  a  very  simple  scaling  law  which  predicts  the  single-scan 
correlation  performance  as  a  function  of  object  density  on  the  focal  plane  and  the  average 
prediction  error  standard  deviation.  These  two  variables  are  time-varying  and  depend  on 
threat/sensor  scenarios.  In  this  section,  we  will  discuss  how  to  obtain  these  two  variables 
from  a  given  set  of  scenario  variables,  together  with  other  related  issues. 

3.1  Cramer-Rao  Type  Bound 

This  sub-section  discusses  how  to  obtain  the  average  prediction  error  standard  devi¬ 
ation  a  defined  by  (2.9).  The  standard  deviation  a  is  determined  by  the  prediction  error 
variance  matrix  5  which  we  are  assuming  (at  least  approximately)  identical  for  every  track- 
measurement  pair.  In  Section  2,  we  have  modeled  the  prediction  error  by  an  additive  gaussian 
random  vector  as  in  (2.2).  For  an  actual  track,  the  measurement  prediction  y  and  the  vari¬ 
ance  matrix  S  are  calculated  from  the  six-dimensional  state  estimate  and  its  error  variance 
prediction  associated  with  the  track.  Unless  all  the  past  measurement  associations  are  given, 
it  is  impossible  to  calculate  or  estimate  the  state  estimation  errors.  To  keep  the  development 
simple,  we  will  first  consider  a  track  with  perfect  past  measurement  association,  which  will 
give  us  a  form  of  upper  bound  on  the  correlation  performance  for  a  given  scan.  The  effects 
of  the  past  data  mis-associations  will  be  discussed  later. 

When  all  the  data  associations  are  known,  estimation  of  the  six-dimensional  state  of  a 
ballistic  object  becomes  a  well-defined  continuous-time  dynamics  discrete-time  measurement 
filtering  problem.  Due  to  the  non-linearity,  there  is  no  analytic  way  of  accurately  predicting 
the  filtering  performance.  Fortunately,  it  is  possible  to  calculate  a  Cramer-Rao  type  lower 
bound  on  the  state  estimation  errors  using  a  method  similar  to  the  one  used  in  [8].  The 
technique  can  be  web  summarized  by  the  following  quotation  from  the  reference  (9): 

“ .  the  extended  Kalman  filter  variance  propagation  equations  linearized 

about  the  true  unknown  trajectory  provides  the  Crame-Rao  lower  bound  to  the 
estimation  error  variance  matrix.” 
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In  other  words,  we  can  solve  the  Riccati-type  differential  equation  using  partial  derivatives 
about  the  true  trajectory  to  obtain  a  lower  bound  on  the  state  estimation  error  variance 
matrix,  from  which  we  can  calculate  a  lower  bound  on  the  measurement  prediction  error 
variances. 

Consider  a  track  which  consists  of  a  set  of  measurements  at  different  times  indexed  as 
{fjk}  (in  general  from  multiple  sensors)  and  assume  that  the  object  being  tracked  is  governed 
by  the  ballistic  dynamical  equation 

-^z(i)  =  f(x(t))  (3.1) 

where  x(t)  is  the  six-dimensional  state  consisting  of  the  3-dimensional  position  and  velocity, 
and  /  is  the  gradient  of  the  potential  due  to  the  earth’s  gravity.  The  track  contains  the 
measurement  at  each  time  tk  as 

yk  =  h(x(tk)1  th)  +  u>k  (3-2) 

with  white  gaussian  noise  sequence  {10*}.  The  function  h(-,  tk)  projects  the  three-dimensional 
object  position  onto  the  sensor’s  focal  plane  at  time 

For  the  sake  of  simplicity,  assume  that  there  is  no  a  prion  information.  Then  the  lower 
bound  Pfc  on  the  state  estimation  (Altering)  error  variance  matrix  at  time  tj,  can  be  written 
as 

Pk1  =  E  (3.3) 

•=1 

where  ¥(•,•)  is  the  fundamental  solution  matrix  satisfying 

j^(t,s)  =  -¥(t,j)^(x(t))  (3.4) 

with  the  boundary  condition  ll,(s,s)  =  /  where  I  is  the  six-dimensional  identity  matrix.  Hk 
in  (3.3)  is  the  matrix  of  the  observation  partial  about  the  true  state  at  time  tk,  i.e., 

Hk  =  ~(x(tk),tk)  (3.5) 

and  Rk  is  the  measurement  error  variance  matrix  for  the  n-th  observation.  Finally  the 
measurement  prediction  error  variance  matrix  for  the  scan  at  time  tk  is  calculated  as 

Sk  =  HhPkHl  +  Rh  (3.6) 


25 


Table  3-1:  Choice  of  State  Space  for  Cramer- Rao  Bound  Calculation 


State  Space 

Dynamics 

Observation 

Cartesian 
Range,  Az,  El 
Orbital  Elem. 

Simple  Nonlinear 
Complicated  Nonlinear 
Linear 

Simple  Nonlinear 
Linear 

Complicated  Nonlinear 

In  the  above  equations,  (3.1)  -  (3.6),  we  have  assumed  that  the  object  state  is  expressed 
in  the  six- dimensional  cartesian  coordinates  with  the  origin  at  the  center  of  the  earth.  How¬ 
ever,  the  results  (3.6)  should  remain  the  same  even  when  we  use  a  different  coordinate 
system.  Table  3-1  shows  the  possible  choices  of  the  coordinate  system  for  the  state  space. 
Appendix  E  describes  the  forms  of  functions,  (/,  h ),  in  each  coordinate  system.  One  of  our 
objectives  is  to  calculate  the  performance  estimates  quickly  given  a  small  set  of  necessary 
threat /sensor  parameters.  The  choice  of  the  coordinate  system  may  be  related  to  this  re¬ 
quirement.  Nonetheless,  we  did  not  find  any  particular  difficulty  in  solving  the  differential 
equations,  (3.1)  and  (3.4),  in  the  standard  six-dimensional  cartesian  coordinate. 

As  a  preliminary  investigation,  we  assume  a  simple  scenario  assuming  a  threat  consist¬ 
ing  of  a  group  of  ballistic  objects  and  a  single  orbiting  sensor.  Figure  3-1  illustrates  the 
threat  and  the  sensor  trajectories.  The  two  circles  in  the  figure  indicate  the  centroid  of  the 
threat  and  the  sensor  position  at  a  given  time.  The  threat  trajectory  is  a  minimum-energy 
trajectory  and  the  sensor  orbit  is  a  low-altitude  circular  one.  This  example  was  constructed 
so  that  the  threat  is  visible  from  the  sensor  for  most  of  the  threat  trajectory.  Figure  3-2 
shows  the  results  of  calculations  of  the  average  prediction  error  standard  deviations  a  as¬ 
suming  perfect  data  association  at  each  scan.  This  figure  shows  the  average  prediction  error 
standard  deviations  calculated  by  (2.9)  and  normalized  by  a  constant  sensor  measurement 
error  standard  deviation  with  three  different  sampling  intervals,  5,  10  and  20  seconds.  Each 
curve  shown  in  the  figure  has  a  peak  around  time  350  sec.  when  the  distance  from  the  sensor 
to  the  threat  group  is  at  a  minimum  so  that  the  prediction  errors  are  magnified,  and  exhibits 
a  pseudo-stationary  behavior  just  before  and  after  this  peak. 

As  seen  in  Figure  3-2,  when  a  track  has  only  a  few  scans  of  data,  the  prediction  errors 
are  very  big,  partially  because  of  the  low  observability  of  the  angle-only  measurements  by  the 
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sensor.  This  poor  prediction  performance  by  tracks  will  result  in  very  poor  correlation  results 
if  a  simple  correlation  algorithm  discussed  in  Section  2  were  used;  in  fact,  correlation  results 
from  such  an  algorithm  may  be  so  bad  as  to  become  meaningless.  Generally,  the  problem 
of  generating  meaningful  tracks  under  such  conditions  is  referred  to  as  the  track  initiation 
problem.  Many  techniques  have  been  developed,  including  the  genera]  multi-hypothesis 
technique  described  in  [l].  Several  techniques  were  developed  specifically  for  ballistic  missile 
defense  systems  (with  which  we  are  concerned  in  this  report),  i.e.,  [10],  [ll]  and  its  refinement. 
From  the  view-point  of  track  purity,  tracks  in  this  track  initiation  phase  may  be  too 
poor  to  be  useful  for  the  user  we  have  in  mind  for  this  report,  i.e.,  an  object  discrimination 
system.  Hence  we  may  consider  a  track  initiation  period  as  a  time  interval  to  produce  a  priori 
information  to  tracks  and  consider  only  subsequent  intervals  for  object  discrimination,  which 
is  more  or  less  consistent  with  the  underlying  assumptions  made  in  our  analysis  described 
in  this  report.  ' 

Under  some  conditions,  however,  it  is  desirable  to  provide  as  many  data  as  possible 
for  object  discrimination  and  even  to  include  the  measurements  in  the  initiation  phase. 
In  such  a  situation,  it  is  reasonable  to  “re-do”  the  data-to-data  correlation  after  tracks 
accumulate  enough  measurements  to  establish  “good”  state  estimates.  One  simple  way  of 
viewing  such  a  process  is  to  consider  track-to-measurement  correlation  using  the  smoothing 
state  estimates  rather  than  the  filtering  state  estimates  as  indicated  by  (3.3).  A  lower  bound 
on  the  estimation  error  variance  for  the  smoothing  state  estimate  can  be  obtained  by  simply 
replacing  the  range  of  summation  in  (3.3)  by  a  larger  interval,  e.g.,  up  to  a  given  K  as 

p;'  =  (3.r) 

»=i 

Correspondingly  we  may  replace  Ph  in  (3.6)  by  P*.  Figure  3-3  shows  an  example  comparing 
the  “filtering”  predictions  with  the  “smoothed”  predictions.  However,  it  is  not  very  clear  if 
the  use  of  the  smoothing  estimation  error  variances  in  this  way  can  be  justifiable  (or  useful) 
or  not.  For  example,  it  is  not  clear  if  we  should  use  the  “smoothed”  state  estimation  error 
variance  to  predict  performance  of  such  “delayed”  correlation  since  the  information  in  the 
measurements  are  already  included  in  the  “smoothed”  state  estimate. 
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Figure  3-3:  Average  Prediction  Error  Standard  Deviation  -  Smoothed  - 
3.2  Object  Density  Estimation 

This  sub-section  discusses  the  calculation  of  the  other  parameter  needed  for  the  simple 
scaling  law,  (2.13)  or  (2.15),  i.e.,  the  object  density  on  the  sensor’s  focal  plane.  We  will  seek 
a  simple  analytic  method  to  calculate  this  quantity  from  a  given  set  of  threat  parameters. 

Consider  N  ballistic  objects  and  let  x»(t)  be  the  6D  state  of  the  i-th  object  at  time  t. 
When  all  the  objects  are  separated  (no  interaction  among  them)  and  the  only  force  applied 
is  the  earth’s  gravitational  force,  all  the  states  Zj(t)  obey  the  same  differential  equation, 
Xi  =  f{xi).  Let  us  define  the  centroid  by 

*(0  =  <3-8> 

•*’  i=i 

and  the  spatial  variance  by 

v(t)  ^  i  £>.-«)  -  «(*))(*(«)  -  mf  (3-9) 

jy  «= i 
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Then,  if  the  distance  ||i<(t)  —  i(t)||  of  each  object  state  from  the  centroid  is  small  enough, 


the  centroid  and  the  variance  follow  the  extrapolation  equations 

=  4  E /(*•(*)) 

1  «•=  1 

~  jj  E  (/(*(*))  +  |^(*(0)(*i(<)  -  *(<))) 

=  /(*(*)) 


(3.10) 


!v(t)  =  |(^g>W*WT-*(')*(<f) 

«  jf  23  (|j(*(t))(Ii(t)  -  i(0)I.(i)T  +  *i(0(*(W  - 
=  ^(*«))V(()+V(()^(*«)T 

(3.11) 

The  spatial  variance  E(t)  on  a  sensor’s  focal  plane  at  time  t  can  then  be  approximated 


as 


s(o  »  (3.12) 

where  h(-,t)  is  the  function  which  maps  a  point  in  the  six-dimensional  state  space  to  a  two- 
dimensional  angular  measurement  at  time  t  on  the  focal  plane  of  a  given  sensor.  The  object 
density  on  the  focal  plane  can  be  calculated  as  a  ratio  of  the  total  number  of  objects  over 
the  area  of  the  elliptic  region  defined  by  the  variance  matrix  S (f),  using  an  appropriate  x* 
value,  as 


/?(<)  = 


N 

7rx2v/det(S(t)) 


(3.13) 


Figure  3-4  shows  an  example  object  density  on  a  sensor’s  focal  plane  as  a  function 
of  time.  The  scenario  used  for  this  calculation  is  the  same  as  the  one  used  in  the  last 
section.  The  initial  spatial  variance  has  been  more  or  less  randomly  chosen.  The  x2  value 
corresponding  to  the  99%-inclusion  probability  was  used  for  (3.13).  The  object  density  in 
the  figure  is  normalized  by  the  sensor’s  measurement  error  standard  deviation.  The  dip  of 
the  curve  in  Figure  3-4  corresponds  to  the  point  where  the  distance  between  the  threat  group 
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Figure  3-4:  Object  Density  on  Sensor’s  Focal  Plane 

and  the  sensor  is  at  a  minimum.  As  seen  in  this  figure,  unlike  the  prediction  error  variance, 
in  general,  there  is  no  apparent  “stationary”  behavior.  The  initial  condition  for  (3.12)  was 
set  to  be 

0  0 

<7„J  o  o 

0  0  crv2  0 

0  0  ov2 

with  a  given  <r„  which  we  can  relate  to  as  the  dispersion  rate  of  a  threat  object  group. 
Needless  to  say,  the  x2  value  in  (3.13)  is  sensitive  to  the  calculation  of  the  object  density  /? 
and  an  appropriate  value  is  yet  to  be  investigated. 

3.3  Effects  of  Mis- Association 

Mis-association  for  a  track  in  a  scan  may  create  significant  measurement  prediction 
errors  in  the  subsequent  scans.  Such  errors  may  be  characterized  as  estimation  bias  in  the 
sense  that  the  estimates  may  be  further  away  from  the  “true”  values  than  being  indicated  by 
the  variance  matrices  accompanying  the  estimates.  Apparently,  the  effects  of  mis-association 
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will  be  propagated  from  scan  to  scan.  They  may  be  magnified  or  be  corrected  in  later  scans. 
For  example,  a  mis-association  in  one  scan  produces  a  bias  which  makes  mis-associations  in 
the  subsequent  scans  more  likely.  On  the  other  hand,  if  the  track  happens  to  be  associated 
with  the  true  measurement  in  the  next  scan,  the  bias  may  be  compensated  and  the  measure¬ 
ment  predictions  may  recover  for  the  subsequent  scans.  Thus  the  mechanism  in  which  the 
effects  of  mis-associations  are  propagated  from  scan  to  scan  is  very  complicated,  and  hence, 
very  hard  to  model.  Recently  it  was  proposed  ((l3j)  to  extend  Markov-chain  modeling  tech¬ 
niques  used  for  estimating  track  initiation  probability  (e.g.,  in  [3],  [4],  [5])  to  model  such  a 
mechanism. 

We  tried  to  incorporate  this  Markov-chain-model  technique  into  our  framework,  but 
no  significant  results  were  obtained.  Difficulties  lies  in  specifying  the  necessary  transition 
probabilities  necessary  for  building  the  Markov-chain  model  to  describe  the  propagation  of 
mis-association  effects.  From  a  limited  number  of  simulations  using  a  Markov-chain  model, 
we  tentatively  concluded  that  the  effects  of  mis-associations  may  not  be  very  significant  in 
a  dense  target  environment.  This  may  be  so  because  even  when  a  mis-association  happens 
the  “wrong  measurement”  must  be  close  enough  to  be  mistaken  so  that  the  effects  may  fall 
within  the  range  of  sensor  measurement  errors.  On  the  other  hand,  if  the  object  density 
as  well  as  false  alarm  density  is  substantially  low,  the  probability  of  having  mis-association 
should  be  very  low.  Thus  there  may  be  an  intermediate  range  of  object  density  where  the 
effects  of  mis-association  may  significantly  affect  the  tracking  performance. 

3.4  Parameterization  of  Threat /Sensor  Geometry 

In  the  calculation  of  the  average  prediction  error  standard  deviation  5  and  the  object 
density  /?,  we  have  utilized  all  the  parameters  defining  the  orbits  of  the  threat  centroid  and 
the  sensor.  In  order  to  have  a  simple  expression,  it  may  be  desirable  to  parameterize  the 
geometric  relationship  between  the  threat  and  the  sensor  orbital  parameters.  In  situations 
where  time  constants  in  threat/sensor  dynamics  are  substantially  longer  than  those  in  our 
environment,  we  may  consider  parameters  such  as  range ,  aspect  angle ,  depression  angle, 
relative  velocity,  etc.  However,  as  seen  in  the  previous  figures  in  this  section,  the  geometric 
relation  may  never  be  stationary.  When  a  particular  sensor  has  a  period  of  blockage,  either 
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due  to  the  earth  field-of-view  blockage  or  the  loss  of  SNR  for  threats  appearing  on  the  earth’s 
surface,  parameterization  becomes  even  more  difficult. 

At  this  point,  it  is  not  clear  at  all  if  we  may  be  able  to  do  any  meaningful  parameteriza¬ 
tion  or  not.  As  indicated  in  Section  3.1.,  use  of  a  special  coordinate  system  for  threat/sensor 
orbital  elements  may  provide  some  insight  on  this  problem. 


4.  TRACK  PURITY  ESTIMATION 


In  previous  sections,  we  have  established  a  method  for  predicting  the  single-scan  data 
association  performance  given  a  set  of  threat /sensor  parameters.  In  this  section,  the  predic¬ 
tion  of  track  purity  based  on  the  single-scan  performance  predictions  will  be  discussed.  We 
will  also  discuss  the  issues  concerning  unresolved  measurements  (CSO’s). 


4.1  Definition  and  Prediction  of  Track  Purity 


As  in  Section  3.1,  let  us  consider  sensor  scans  (in  general  from  multiple  sensors)  indexed 
by  time  tj,.  Suppose  that  we  are  given  a  set  K  of  scan  indices  and  that  we  wish  to  define 
the  purity  of  a  given  track.  For  example,  A  is  a  set  of  scan  indices  for  several  consecutive 
scans  from  one  of  the  sensors  used  in  a  tracking  system.  It  is  natural  to  define  the  purity  as 


^  ,  -  .  No.  of  scans  with  correct  measurement  assignments 

Track  Punty  =  - — - - - 

Total  no.  of  scans  in  K 


(4.1) 


An  immediate  problem  with  this  definition  is  the  treatment  of  unresolved  (merged)  mea¬ 
surements  and  mis-detections.  The  problem  of  limited  sensor  resolution  will  be  discussed 
in  the  next  two  sections.  On  the  other  hand,  it  is  not  quite  clear  even  how  to  treat  failed 
detections.  If  an  object  is  not  detected  but  associated  with  a  measurement,  it  is  clearly  a 
mis-association.  But  if  it  is  correctly  recognized  (as  undetected),  it  is  not  clear  if  we  should 
count  it  as  a  “correct  association”  or  not.  In  this  report,  by  assuming  that  the  detection 
probability  is  high  enough,  it  was  decided  not  to  include  mis-detections  in  the  discussion  of 
track  purity. 

The  purpose  of  this  section  is  to  produce  a  reasonable  prediction  on  the  average  track 
purity  given  threat/sensor  parameters.  For  this  purpose,  consider  &  fictitious  binary  stochas¬ 
tic  process  {m*}  with  discrete  index  k  in  the  set  K.  The  event,  m*  =  1,  represents  a 
track-measurement  mis-association,  i.e., 


Prob.{mfc  =  0}  =  Pc(j3k,ok)  =  exp(-7r  ^ko\) 


(4.2) 


where  Pc{0h,<?k)  is  the  (average)  probability  of  a  track  being  correctly  associated  with  a 
measurement  at  the  scan  k  as  a  function  of  the  object  density  and  the  average  prediction 
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error  standard  variation  o/,.  In  fact,  the  function  Pc  is  the  scaling  law  described  by  (2.13) 
in  Section  2.2.  The  total  number  of  mis-associations  in  a  given  set  of  scans  is  then  simply 

M[ff]  =  £  m*  (4.3) 

keK 

We  may  view  the  track  purity  defined  by  (4.1)  for  each  track  as  a  sample  of  the  random 
variable1,  M [K]/#(K),  and  hence,  calculate  the  average  track  purity  as 

m  rc(A,»*)  (4.4) 

keK 

When  max*{l  —  Pc(flk,<Tk)}  is  small  enough,  we  may  apply  the  law  of  small  numbers  (See, 
e.g.,  [14])  and  show  that  we  may  approximate  the  distribution  of  the  random  variable  M[JT] 
by  a  Poisson  distribution  with  the  mean 

I'M  =  S  C1  ~  pc(0k°k))  »  (4-5) 

fceff 

4.2  Treatment  and  Prediction  of  Unresolved  Measurements  (CSO’s) 

The  first  question  we  must  answer  concerning  unresolved  measurements  (CSO’s)  is  how 
we  should  count  such  measurements  in  calculating  the  track  purity.  When  we  view  tracking 
as  a  state  estimation  problem  with  multiple  objects,  we  may  consider  a  correlation  result  as 
correct  if  a  track  is  associated  with  an  unresolved  measurement  with  multiple  origins  one  of 
which  is  the  “true”  origin  of  the  track.  In  general,  an  unresolved  measurement  contains  a 
bias  when  it  is  considered  as  a  measurement  for  a  given  track,  thereby  somewhat  discounting 
the  “correctness”  of  the  association. 

As  mentioned  before,  the  primary  concern  of  this  report  is  the  track  purity  seen  from 
an  object  discrimination  system.  When  two  or  more  objects  create  a  single  detection,  it  is 
generally  extremely  difficult  to  separate  out  the  contribution  for  each  separate  object  from  the 
merged  detection.  At  least,  this  is  not  clear  when  examining  a  functional  model  developed  for 
simulation  purposes  (  [15]).  If  indeed  it  is  not  possible  to  extract  useful  information  for  object 
discrimination  from  unresolved  measurements,  then  it  is  appropriate  to  count  the  unresolved 
measurements  as  impure  measurements  which  are  as  bad  as  mistaken  measurements.  When 
‘By  #(4),  we  mean  the  cardinality,  i.e.,  the  number  of  elements,  of  a  set  A. 
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Figure  4-1:  Probability  of  Measurement  Merging  (1) 

we  taking  this  point  of  view,  we  should  simply  discount  the  track  purity  by  the  probability 
of  the  detection  (corresponding  to  a  track  in  question)  being  unresolved. 

We  will  base  our  calculation  of  the  probability  of  measurements  being  unresolved  on  the 
functional  model  described  in  [15].  Consider  two  objects  represented  by  two  two-dimensional 
vectors,  yi  and  y2.  According  to  [15],  the  probability  of  measurements  from  these  two  objects 
being  unresolved  can  be  modeled  by  a  function  pm  of  the  distance  d  =  ||yj  —  y2||  between 
the  two  objects.  A  good  model  for  pm  is 

1  if  0<d<6/2 

Pm(d)  =  *  2  — - —  else  if  6/ 2  <  d  <  6  (4.6) 

b 

0  otherwise 

which  is  shown  in  Figure  4-1,  where  6  is  the  sensor  resolution.  With  this  model  and  a  set  of 
assumptions,  we  can  show  (See  Appendix  F.)  that  an  appropriate  expression  of  track  purity 
considering  unresolved  measurements  is 

■jrhn  £  -  P"<M) 

)  keK 


(4.7) 


Figure  4-2:  Probability  of  Measurement  Merging  (2) 
where  Pm  (A,)  is  the  probability  of  a  track  merged  with  another  track  and  is  expressed  as 
Pm  (At)  =  1  -  ^erf  (^J 2^6 )  -  erf  ^2x^0  j  (4.8) 

with  5  being  the  sensor  resolution.  (See  Appendix  F  for  the  derivation.) 

Figure  4-2  shows  the  measurement  merging  probability  Pm  as  a  function  of  object 
density.  The  curve  in  the  figure  is  the  plot  of  (4.8)  while  each  point  (small  filled  circle)  is 
obtained  by  500-run  Monte  Carlo  simulations  with  100  simulated  objects  for  each  run.  In 
the  figure  the  object  density  is  normalized  by  the  sensor  resolution  cell  area,  i.e.,  62.  As  seen 
in  this  figure,  the  expression  (4.8)  provides  us  with  an  accurate  estimate  on  the  probability. 
Figure  4-3  shows  the  probability  of  correct  association  using  different  object  density  values 
and  the  deterioration  of  the  performance  due  to  the  unresolved  measurements.  In  this  figure, 
one  curve  (broken)  represents  the  expression  ignoring  the  CSO’s,  i.e.,  (4.4),  and  the  other 
(solid)  the  one  considering  them,  i.e.,  (4.7).  Each  point  in  the  figure  is  the  actual  correlation 
performance  obtained  by  simulations  with  unresolved  measurements.  In  the  simulations,  a 
functional  model  described  in  [15]  was  used.  This  functional  model  takes  the  deterioration 
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of  measurement  accuracy  due  to  neighboring  measurements  into  account.  As  seen  in  this 
figure,  (4.7)  accurately  predicts  the  correlation  performance  with  unresolved  measurements. 

4.3  Simulation  Results 


The  test  scenario  of  Figure  3-1,  consisting  of  a  threat  object  group  and  an  orbiting 
sensor,  was  used  to  compare  the  analytic  prediction  of  the  average  track  purity  with  simulated 
tracking  results.  A  simple  tracking  algorithm  was  used  for  this  purpose.  In  particular,  track 
initiation  problems  were  avoided  by  using  relatively  “good”  a  priori  information. 

The  tracking  algorithm  may  be  summarized  as  follows: 


1.  For  each  scan,  first  the  object  state  distributions  of  all  the  existing  tracks  are  extrap¬ 
olated  to  the  time  of  the  scan,  and  then  measurement  predictions  are  made  for  each 
track  using  the  sensor  position  at  the  time. 


2.  Then,  for  each  existing  track,  an  association  list  consisting  with  potential  measurements 
and  the  track-measurement  association  x 2  values  is  created. 
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3.  Munkres  algorithm  is  then  used  to  obtain  the  optimal  assignment  in  terms  of  track- 
measurement  pairs. 

4.  The  tracks  which  are  not  assigned  to  any  measurement  are  then  examined.  For  each 
of  such  tracks,  if  a  measurement  is  found  such  that  the  %2  value  is  within  a  given 
threshold,  the  track-measurement  pair  is  added  to  the  assignment  list  made  in  the 
above  step.  (This  step  is  to  account  for  tracks  which  have  been  resolved  in  the  past 
scans  but  not  resolved  in  the  current  scan.) 

5.  The  object  state  distributions  of  all  the  tracks  included  in  the  assignment  list  are  then 
updated  by  the  assigned  measurements. 

6.  All  the  measurements  which  are  not  assigned  to  any  track  are  then  examined.  For  each 
of  such  measurements,  a  track  for  which  the  x2  value  is  minimized  is  then  searched. 
If  the  minimal  x2  value  is  within  a  given  threshold  value,  then  the  track  is  split  using 
this  measurement. 

7.  The  remaining  measurements  are  then  examined  as  potential  sources  to  initiate  tracks 
together  with  the  a  priori  object  state  distribution  attributed  to  the  threat  group  with 
a  simple  thresholding  scheme. 

8.  Then  all  the  tracks  are  checked  and,  according  to  a  simple  rule,  redundant  or  “bad” 
tracks  are  eliminated. 

We  have  conducted  two  simulations;  one  with  all  the  measurements  being  resolved  and 
the  other  one  with  unresolved  measurements,  to  clarify  the  effects  of  CSO’s.  Unfortunately, 
due  to  time  and  other  resource  constraints,  only  a  single  sample  run  for  each  case  was 
conducted,  and  hence,  we  cannot  conclude  the  validity  of  our  analytic  models  beyond  sample 
variations. 

The  parameters  used  in  the  scenario  are  the  same  as  the  ones  described  in  Section 
3.1.  Figure  4-4  shows  the  estimated  average  track  purity  calculated  according  to  (4.4)  (i.e., 
without  considering  unresolved  measurements)  for  three  different  sampling  intervals,  5,  10 
and  20  seconds.  Figure  4-5  shows  the  Cramer- Rao  type  lower  bound  on  the  prediction  error 
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Figure  4-4:  Estimated  Average  1>ack  Purity  As  a  Function  of  Time  (1) 


Figure  4-5:  Prediction  Error  Standard  Deviation  (Simulated  Data) 
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Figure  4-6:  Average  Track  Purity  (Estimated  vs  Simulated)  (1) 


standard  deviations  (used  to  calculate  the  curves  in  Figure  4-4).  Figure  4-5  also  plots  the 
average  prediction  error  standard  deviations  calculated  from  actual  tracks  in  the  simulation. 
The  “lower  bounds”  in  Figure  4-5  are  greater  than  the  actual  standard  deviations.  This 
is  so  because  the  lower  bounds  were  calculated  assuming  that  no  a  priori  information  is 
available.  Thus  the  differences  between  the  analytic  curves  and  the  actual  averages  (from  the 
simulation)  show  the  effects  of  the  a  priori  information  we  have  assumed  in  the  simulation. 
Figure  4-6  shows  the  comparison  of  the  average  track  purity  predicted  by  our  analysis  with 
the  one  obtained  from  the  simulation,  with  three  different  sampling  intervals.  Some  persistent 
difference  between  the  analytic  results  and  simulation  results  is  seen  in  this  figure.  As 
mentioned  before,  however,  the  simulation  results  are  only  from  one  sample  run,  and  hence, 
it  may  not  be  possible  to  draw  any  definite  conclusion  from  these  results,  except  for  the  fact 
that  the  simulation  results  are  generally  close  to  our  analytically  predicted  curves. 

In  the  simulations  (using  three  different  sampling  intervals  and  the  results  of  which 
are  shown  in  Figures  4-4  to  4-6),  all  the  measurements  are  resolved  and,  assuming  a  con¬ 
stant  SNR,  a  fixed  sensor  measurement  error  variance  was  assumed.  Another  simulation 
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Probability  of  Measurement  Merging 


\ 

Figure  4-7:  Probability  of  Measurement  Merging  As  a  Function  of  Time 


was  conducted  with  unresolved  measurements.  The  functional  model  described  in  [15]  was 
used  for  sensor  measurement  simulation  assuming  a  constant  SNR.  First,  the  probability 
of  measurement  merging  was  calculated  as  a  function  of  time,  using  (4.8).  The  results  are 
shown  in  Figure  4-7.  A  dip  of  the  curve  in  the  figure  corresponds  to  the  time  when  the 
distance  between  the  threat  group  and  the  sensor  is  minimum  and  so  is  the  object  density. 
Figure  4-8  shows  two  average  track  purity  curves  as  functions  of  time,  which  are  predicted 
tom  our  analysis;  one  with  unresolved  measurements  and  the  other  assuming  all  the  mea¬ 
surements  are  resolved.  The  difference  between  the  two  curves  is  due  to  the  probability  of 
measurement  merging  shown  in  Figure  4-7.  Figure  4-9  compares  the  track  purity  predicted 
(4.7)  (considering  the  unresolved  measurements  as  “impure  measurements”)  with  the  results 
from  a  simulation  run.  Although  there  is  some  initial  disagreement  between  the  two  results, 
in  general,  the  predicted  track  purity  is  very  closed  to  that  obtained  from  the  simulation, 
ESPECIALLY  as  the  number  of  scans  is  accumulated,  .  When  evaluating  the  simulation 
results,  all  the  past  measurement  assignments  of  each  track  at  various  stage  are  examined 
and  the  majority  of  the  origins  of  the  past  measurements  is  then  used  for  the  “true"  origin  of 


42 


Figure  4-8:  Estimated  Average  'Rack  Purity  As  a  Function  of  Time  (2) 


Figure  4-9:  Average  Track  Purity  (Estimated  vs  Simulated)  (1) 


•  ,  ,  as  the  retio  of  the  number  of  measurements  from 

«*  ^  ^ ;  ;PMt *  *.  ^ 
the  “  W  origin  »»«  *h«  ‘°tal  numb”  P 


5.  CONCLUSIONS 


A  very  simple  analytic  model  was  developed  for  predicting  the  tracking  performance 
in  terms  of  average  track  purity,  as  a  function  of  threat/sensor  parameters.  This  model 
is  based  on  a  scaling  law  which  describes  the  single-scan  track-to-measurement  correlation 
performance  as  a  function  of  the  object  density  on  the  sensor’s  focal  plane  and  the  average 
prediction  error  standard  deviation.  Fairly  extensive  Monte  Carlo  simulations  were  con¬ 
ducted  to  validate  the  analytic  expression  of  the  single-scan  correlation  performance.  It  was 
found  that  the  Monte  Carlo  simulations  support  the  predictions  by  our  analysis  very  well 
in  a  reasonable  range.  Assuming  a  very  simple  tracking  algorithm  and  relatively  accurate  a 
priori  information,  a  limited  number  of  small-scale  tracking  simulations  were  also  conducted 
to  verify  the  ability  of  our  analytic  models  to  predict  track  purity. 

5.1  Achievements 

The  achievements  of  our  research  efforts  can  be  summarized  as  follows: 

1.  Simple  expressions  were  developed  to  assess  single-scan  track-to-measurement  correla¬ 
tion  performance  and  to  estimate  the  performance  of  a  given  set  of  sensors  against  a 
given  threat,  in  terms  of  average  track  purity. 

2.  Fairly  extensive  Monte-Carlo  simulations  support  the  simple  analytic  model  for  pre¬ 
dicting  single-scan  track-to-measurement  correlation  performance. 

3.  The  effects  of  eccentricity  of  the  prediction  error  ellipses  were  examined  to  verify  the 
calculation  of  correct  association  probability. 

4.  Uneven  correlation  performance  due  to  A  non-uniform  object  density  was  investigated. 

5.  The  effects  of  false  alarms  were  investigated.  A  simple  extension  to  the  analytic  model 
was  made  to  account  for  such  effects.  Fairly  extensive  Monte-Carlo  simulations  has 
shown  the  validity  of  this  extension. 
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6.  Cramer- Rao  type  lower  bounds  on  state  estimation  error  variance  matrices  as  functions 
of  a  threat  centroid  trajectory  and  a  sensor  orbit  were  calculated  to  predict  the  track 
accuracy  through  tracking  processes. 

7.  A  simple  model  to  describe  object  density  on  the  sensor’s  focal  plane  was  derived. 

8.  The  average  track  purity  as  a  function  of  time  was  calculated. 

9.  The  effects  of  unresolved  measurements  were  analyzed.  In  particular,  a  simple  ex¬ 
pression  for  the  probability  of  measurement  merging  in  terms  of  sensor  resolutions 
and  object  densities  was  developed  and  incorporated  into  the  calculation  of  the  track 
purity. 

10.  Small-scale  Monte-Carlo  tracking  simulations  have  demonstrated  the  validity  of  apply¬ 
ing  the  simple  analytic  model  to  actual  tracking  system  performance  analysis. 

5.2  Topics  of  Future  Research 

The  issues  deferred  to  future  research  are  summarized  as  follows: 

1.  Further  development  of  analytic  performance  models 

a.  Refine  performance  models:  As  far  as  single-scan  models  are  concerned,  simple 
analytic  models  developed  in  the  period  covered  by  this  report  have  been  proved 
to  be  capable  of  accurately  predicting  correlation  performance.  There  is,  however, 
room  for  refinement  in  modeling  multi-scan  performance.  In  particular,  the  effects 
of  mis-associations  on  the  correlation  performance  in  the  subsequent  scans  should 
be  investigated.  The  effects  of  unresolved  measurements,  mis-detections  and  false 
alarms  need  to  be  addressed  in  greater  details. 

b.  Develop  correlation  stochastic  models:  From  the  point  of  view  from  an  object 
discrimination  system,  correlation  statistics  of  track  purity  (both  as  a  time  series 
and  as  a  sample  set  at  one  time)  may  be  as  important  as  the  average  track  pu¬ 
rity.  Such  correlation  statistics  may  be  obtained  by  analyzing  interaction  between 
correlation  performance  and  prediction  accuracy  by  tracks. 
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2.  Validation  of  analytic  performance  models: 

a.  Validate  by  more  extensive  simulations:  Only  a  small  number  of  test  scenarios 
were  conducted  for  the  study  described  in  this  report.  More  studies  extensive 
multi-scan  (possibly  multiple-sensor)  simulations  must  be  conducted.  Through  such 
studies  key  parameters  in  threat /sensor  scenarios  should  be  identified.  Performance 
analysis  should  be  extended  to  cover  more  difficult  tracking  situations  with  required 
accuracy,  such  as  performance  in  autonomous  track  initiation  phases. 

b.  Refine  analytic  model  using  simulation  results:  Through  such  simulations, 
we  may  find  needs  for  further  refinement  of  analytic  models  or  even  for  development 
of  yet  other  models. 

3.  System  analysis: 

a.  Develop  models  to  parameterize  sensor/threat  geometry:  Since  the  number 
of  parameters  which  completely  describe  a  given  threat/sensor  scenario  is  not  small 
and  interaction  among  them  is  not  simple,  there  is  a  need  to  identify  key  parameters 
which  dominate  the  tracking  performance. 

b.  Develop  performance  trade-off  curves:  For  system  design  purposes,  perfor¬ 
mance  trade-off  curves  may  be  found  very  useful.  Analytic  models  developed  in 
the  current  phase  or  the  next  phase  of  research  may  be  used  to  produce  several 
trade-off  curves  useful  for  sensor  system  design. 

c.  Investigate  effects  of  track-initiation:  The  impact  of  track  initiation  perfor¬ 
mance  or  of  hand-over  information  should  be  further  investigated.  In  particular, 
the  effects  of  hand-over  information  or  the  lack  of  it  are  very  important  in  designing 
the  sensor  systems  for  ballistic  missile  defense. 
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A.  Single-Scan  Correlation  Problem 


This  appendix  shows  how  a  correlation  problem,  in  particular,  a  track-to-measurement 
correlation  problem,  can  be  formulated  as  a  classical  assignment  problem.  This  fact  has 
been  treated  as  a  trivial  matter  in  the  past  in  a  majority  of  the  multiple  target  tracking 
literature  and  seldom  justified  analytically.  However,  it  might  be  worthwhile  to  examine  the 
underlying  assumptions  which  support  this  “fact.” 

As  in  Section  2.1,  given  a  positive  integer  N,  we  assume  three  systems,  Y  =  (j/i)£Ln 
Y  —  and  Z  =  (z,)^,  of  random  vectors,  and  a  random  permutation  ir  on  the  set 

{l, ...,  N},  such  that 

P(Y,r,Z,i r)  =  P(Y\Z,x)P(Z\Y)P(ir)P(Y)  (A.l) 


with  P(- )  and  P(-|-)  being  the  usual  shorthand  notation  for  the  a  probability  (unconditional 
and  conditional)  distribution  or  its  density.  In  other  words,  we  assume  a  Markov  chain, 
{y  — ♦  Z  — *  y},  and  that  tt  is  independent  of  Y  and  Z.  The  transition  probabilities  are 
given  as 

P  fw",  6  n  d-‘  I  (*>l)  =  n c(*  - ii ;  (A.2) 

\  i=l  )  i=\  ' 


and 


p  f(y»)£i  e  n  dv<  |  (z<)ili>7r)  =  II  G(y*(o  -  z< ;  R*(i))dyi  (A-3) 

v  »=i  /  »=i 

where  (Qi)^  and  (Rj)jt j  are  the  two  jV-tuples  of  non-negative  definite  symmetric  matrices, 
and  G( •;  V)  is  the  probability  density  function  for  the  zero-mean  gaussian  random  vector  with 
the  variance  matrix  V . 

The  problem  can  then  be  described  as  a  problem  for  determining  an  optimal  estimate 
of  the  random  permutation  7r  given  available  data,  (Y,  Y).  The  maximum  likelihood  (ML) 
estimate,  if,  is  an  estimate  which  maximizes  the  likelihood, 


P{Y,Y\t)  =  P{Y\Y  ,it)P{Y) 

=  /  P{Y\Z,n)P{Z\Y)dZP(Y) 


n 


(A.4) 


=  ^C(*)£L.)  n  -  * ;  *-<o) 


i=l 


with  Sij  =  Qi  -f  Rj.  Therefore,  given  (Y,Y),  a  maximum  likelihood  estimate  if  should 
minimize  the  correlation  cost  function, 


N 


A *)  =  J2  !ly-(.)  -  y.lls- 


i  =  l 


(•) 


(A.5) 


assuming  that  all  S^,’ s  are  positive  definite.  Moreover,  if  we  assume  that  all  the  realizations 
of  the  random  permutation  7r  are  equally  likely,  the  maximum  likelihood  estimate  is  also  a 
maximum  a  posteriori  probability  (MAP)  estimate,  since 


P(Y,Y) 


P{Y,  Y)  N\ 


(A.6) 
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__  We  should  note  that  this  conclusion  is  totally  independent  of  the  a  priori  distribution, 
P(Y).  An  alternative  model  can  be  constructed  using 

P(Y,Y,Z,tc)  =  P(Y\Z,n)P(Y\Z)P(z)P(Z)  (A.7) 

instead  of  (A.l),  and 

p  ((yi)iLi  £  II  din  =  II  -  2«  5  QMh  (a. 8) 

\  t=i  /  t=i 

instead  of  (A.2).  With  this  model,  we  have 

P(Y,Y\*)  =  /  P(Y\Z,*)P(Y\Z)P(Z)dZ  (A. 9) 

in  place  of  (A.4).  As  we  see  in  (A. 9),  if  we  can  consider  P(Z)  as  constant,  e.g.,  when 
Z  =  (zOili  is  i-i.d.  with  a  common  uniform  distribution  with  a  large  enough  support,  the 
maximization  of  (A. 9)  is  reduced  to  the  minimization  of  the  association  function  (A. 5). 

Thus  either  of  these  two  different  models  actually  support  the  association  cost  function 
(A.5).  The  difference  between  these  two  models  can  best  described  by  the  following  two 
assumptions:  \ 

[.4.1]  are  independent  from  each  other. 

(which  is  stated  in  Section  2.1.),  and 

[42]  m1,...,Tn/srj  are  independent  from  each  other. 

By  a  casual  observation,  the  assumption  [42]  might  look  more  “natural”  than  [4l].  As  seen 
in  the  next  appendix,  Appendix  B,  however,  the  first  assumption  will  make  further  analysis 
much  easier  than  the  second  one.  However,  we  should  remember  that  the  model  described 
in  Section  2.1  and  this  appendix  was  to  model  a  track-to-measurement  association  process 
in  a  multi-target  tracking  system.  We  will  see  in  the  following  that  with  such  a  background 
the  first  assumption  appears  more  appropriate. 

Let  us  consider  a  more  or  less  idealized  situation  in  a  multi-target  tracking  system.  Sup¬ 
pose  that  we  have  been  tracking  N  objects  perfectly,  i.e.,  that  the  N  targets  are  independent, 
each  object  i  has  been  uniquely  associated  with  object- wise  independent  information  Ii,  and 
we  have  exactly  N  tracks  so  that  the  t-th  track  is  the  object  state  distribution  conditioned 
by  Ii.  Let  z,  be  the  “true”  projection  of  the  t-th  object  state  onto  the  measurement  space 
and  let  us  assume  that 

P(zijli)  =  G(zi  -  yi  ;  Qi)  (A. 10) 

which  means1 

y,  =  £(zi\Ii)  and  Q,  =  VfclJi)  (A.ll) 

Assuming  the  independent  permutation  7r,  the  a  posteriori  probability  of  the  permutation 
becomes 

P(*\I,Y)'=  (A-12) 

’By  £(T)  and  V(  j-),  we  mean  the  conditional  mean  and  variance  operators. 
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where  I  =  Assuming  that  the  random  permutation  ir  is  totally  random,  in  the  sense 

of  either  the  MAP  or  the  ML,  the  track-to-measurement  association  problem  becomes  an 
optimization  problem  to  maximize  the  likelihood  function  P(Y\x,I).  The  assumption  of 
independent  measurement  errors  together  with  (A. 11)  leads  to 

N 

P{Y  I  *,/)  =  n  G  (»«(.)  -  Vi  5  $.(.•))  (A. 13) 

«=i 

Apparently,  the  maximization  of  this  quantity  is  equivalent  to  the  minimization  of  the  cor¬ 
relation  association  cost  function  (A.5). 

Since  y,  is  the  conditional  expectation,  it  is  orthogonal  to  the  “prediction”  error,  z;  —  yt-, 
in  the  sense, 

£  ((*.-  -  yi)TVi  |  /»)  =  0  (A. 14) 

We  may  stretch  this  fact  further  to  assume  that  the  measurement  prediction  y,  is  independent 
to  the  prediction  error 

mi  =  Zi  -  Hi  (A. 15) 

which  should  be  zero-mean,  gaussian  with  the  variance  matrix  Qi.  This  observation  seems 
to  support  the  assumption  [Xl]  more  than  the  other  one  [.42]. 


B.  Calculation  of  Correct  Association  Probability 

This  appendix  describes  detailed  derivations  of  the  calculation  of  the  correct  association 
probability  described  in  Section  2,  and  in  particular,  derivations  of  eqns.  (2.7),  (2.10)  with 
(2.9)  and  the  main  result  of  Section  2,  i.e.,  eqn.  (2.13). 

As  in  Section  2,  we  will  first  consider  the  cases  where  N  =  2  and  the  random  transpo¬ 
sition  (i.e.,  mis-association).  Let  us  denote  the  event  in  which  the  two  tracks  are  associated 
with  correct  measurements  by  E\2.  Then  this  event  can  be  identified  with  the  following 
condition: 

A J  IS'  J({(1,1),(2,2)})-J({(1,2),(2,1)}) 

=  llvi  ~  Sills-.  +  IlSi  -  jillj-i  -  llvi  -  Mis-'  -  ilVl  -  Siiis- 

y  33  2  31  2  12  /p 

=  II”!  -  ™i||(S-i-S-‘)  +  lln2  -  m2|l(5-l-S-1)  "  lift  ~  y*ks£+s£) 

+  2 (5f21(n2  -  m2)  -  S22(n 2  -  m!))r(yi  -  y2) 

<  0 

We  may  calculate  the  (conditional)  characteristic  function,  V’(C)  =f  £(ejAJf  lyi ,  3/2),  of  the 
random  variable  AJ  as  ' 

V’(C)  =  (det(27r5n)det(27r522))“1/2  •  exp  (— j||yi  -  ft||(s-‘+s-‘))  ’ 

Jnt  exP  {3  (llui||(s-‘-5-1)  +  llu2||(sJ--J‘-s-‘)  +  2(S12U2  +  S^Ui)T(yi  -  2/2)))  (B.2) 

-exp  (~i  (IM#  +  Mss*))  du'du* 

which  can  be  expressed  in  a  closed  form.  Then,  using  some  form  of  numerical  integration, 
possibly  by  an  FFT  program,  we  can  calculate 

Prob.{^2|yl,ya}  =  Prob.{AJ  <  0|yj,y2}  =  /°  /°°  e-*1>{Qd(dz  (B.3) 

Z7T  J —00  ^  —00 

The  results  of  such  a  calculation  is  described  in  [6].  Although  one  can  make  several  interesting 
observations  from  results  of  such  numerical  analysis,  it  is  not  easy  to  extend  the  results 
beyond  the  simplified  two-object  cases  or  to  relate  the  results  to  key  parameters  in  track-to- 
measurement  association  problems. 

On  the  other  hand,  if  we  assume 

S  d^f  S„  =  S22  =  Sl2  =  S21  (B.4) 

(B.l)  becomes 

AJ  = 


< 


— 2(3/1  -  y2)TS  '{yi  -y2) 

2(771!  -  m2  -  (na  -  n2)  -  {{/!  -  y2))T S-1^  -  y2) 
2 {10  -  || yi  -y2||j-i) 

0 
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where  w  =  (mi  —  m2  —  n2  +  n2)rS-1(yi  ~  yz)-  Since  we  have  assumed  that  m^’s  and  n^’s  are 
independent  of  in's,  w  is  gaussian  given  yi  and  y2  with 


and 


£(H*i»!ft)  =  0 

V(w\yuy2)  =  2||yi  —  3/2III-1 


(B.6) 

(B.7) 


from  which  eqn.  (2.7)  follows. 

In  order  to  calculate  the  unconditional  probability,  Prob.(£'i2),  we  need  a  mathematical 
model  for  (yi,y2).  For  this,  let  us  assume  that  the  two  random  vectors  are  independent  and 
have  a  common  distribution  which  is  uniform  on  a  disk  with  the  center  at  the  origin  and  with 
a  given  radius  r.  Before  considering  Prob.(£?i2),  let  us  consider  a  more  restricted  problem 
of  finding  a  conditional  probability  Prob.(f?u|yi  =  0).  Since  y2  is  independent  of  y\ ,  (even 
with  the  condition  y2  =  0)  y2  is  uniformly  distributed  on  the  disk.  Thus  it  follows  from  (2.7) 
that 

Prob(£12  |  Hi  =  0)  =  (at2)-1  /  erf  (P3^s  1 

■/{||vj|j<’-}  V 


y/2 


dy  2 


=  1  -  (it**2)-1  I  erfc  (  ^2|j|  -~\  dy2 

Aiifeii<->  \  V2  J 


(B.8) 


where  erfc(£)  =  1  -  erf(£).  Without  loss  of  generality,  we  can  assume  S  =  Diag.[a2, 62], 
with  a  >  b.  If  5  is  not  diagonal,  we  can  always  diagonalize  it  by  properly  rotating 
the  coordinate  without  changing  the  integral.  Thus,  with  a  coordinate  transformation, 
y2  =  [p  cos  0,  p  sin  0]T ,  we  have 


Jm 


erfc 


1ft  IIS-' 

V2 


-)  dy2  =  ciic(U(0)p)d0pdp 

=  Jo  llo  erfc(^(^)F)^ 

=  [ 3*U(0 )"2  I{U{0)r)  d0 
Jo 


(B.9) 


where 


and 


U{9)  = 


1  II* 


2II5-1 


y/2  ||fc||  \ 


1  /  cos2 

2  [~a> 


20  sin2  O' 
+  ~, 


(B.10) 


I(x)  *  j\dc(0(d( 

[(x*  —  l)erfc(x) 


1  1 
=  4  +  2 


y/2w 


xexp 


(-H 


(B.n) 


Figure  B-l  show  the  function  /(•)  defined  by  (B.ll).  As  we  see  in  this  figure,  this 


I(X) 


Figure  B-l:  Integral  J(i) 


function  converges  to  its  asymptotic  value,  1/4,  approximately  at  3  (corresponding  to  the 
3 -a  point).  On  the  other  hand,  we  have 

fin  r*/2  d6 

/  U{e)-2dO  =  8  /  — - --f--  =  4tt ab  (B.12) 

Jo  Jo  a-2  cos2(0)  +  f>_2  sin2(0) 

Since  y/2U{6)  >  1  fa  for  all  0,  it  follows  from  eqns.  (B.8)  through  (B.12)  that 

Prob.(JE?n  |  yx  =  0)  =  1  -  ^  +  o(r-3e-(1/4)(r/o)’)  (B.13) 

Now  let  us  consider  the  unconditional  probability:  It  follows  from  (2.7)  and  the  i.i.d.- 
uniform  assumption  that 


Prob.(£12)  =  (x r2)~2  f  [  erf  dyx 

K  V  '  h\\ii\\<r)hm\<r)  V  \/2  / 


(B.14) 


where 

Using  the  same  argument  as  the  one  which  has  lead  us  to  (B.13),  we  may  conclude  the 
function  W(yx)  can  be  well  approximated  by  a  constant  function,  i.e.,  JP(0)  which  is  nothing 
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but  (B.13),  if  Hjiiil  <  r  or  practically  \\yi\\  <  r  —  3a.  In  other  words,  except  for  the  cases 
where  yi  locates  on  the  peripheral  of  the  disk ,  we  can  approximate  W(yi)  =  Prob^E^lyj) 
by  W{ 0).  Since  the  area  of  the  periphery  is  proportional  to  r  while  the  whole  area  of  the 
disk  is  proportional  to  r2,  when  averaging  over  yx,  we  can  ignore  this  edge  effect.  Noting 
det(S)  =  a2b2,  we  have  proved  that,  for  large  enough  r,  we  have  (2.10)  with  (2.9). 

Now  suppose  that  the  number  N  of  objects  is  a  random  integer  having  a  Poisson  dis¬ 
tribution  with  mean  v  and  assume  that  the  approximation 

is  valid  where  Ee  is  the  event  where  a  track  is  associated  with  the  correct  measurement 
(following  the  reasoning  described  in  Section  2.1).  Then,  defining  Ee  such  that 
Prob.(2?c|./V’  =  0)  =  Piob.(Ec\N  =  1)  =  1,  we  have 


Prob.(i?c  |  N  =  n) 


Prob.(£c)  « 


(B.17) 


which  is  actually  eqn.  (2.13).  The  last  step  of  the  above  derivation  uses  the  approximation 
based  on  the  assumption  j<r.  The  is  the  object  density  defined  by  (2.14). 


C.  Effects  of  Non-Uniform  Local  Object  Density 


In  Appendix  B,  we  derived  a  simple  scaling  law  to  determine  the  probability  of  correct 
association  as  a  function  of  object  density  and  track  measurement  prediction  errors.  To  do 
so,  we  have  assumed  objects  are  distributed  uniformly  on  a  given  region  of  the  sensor’s  focal 
plane.  This  uniformity  assumption  may  not  validly  reflect  reality.  For  example,  if  objects  are 
distributed  uniformly  within  a  ball  in  a  three-dimensional  space,  apparently  their  projections 
onto  a  (2-dim.)  plane  will  not  have  uniform  density. 

In  this  appendix,  we  will  calculate  the  probability  of  correct  association  for  a  group  of 
targets  distributed  spatially  with  a  gaussian  distribution  instead  of  uniform  distribution.  As 
we  did  in  Appendix  B,  we  will  first  consider  the  transposition  in  the  two-object  model.  We 
assume  that  the  two-dimensional  gaussian  vectors,  yi  and  1/2,  are  independent  but  share  a 
common  zero  mean  gaussian  distribution  with  a  variance  matrix  V.  In  order  to  examine 
the  effect  of  this  “uneven”  (or  non-uniform)  spatial  distribution,  we  would  like  to  calculate 
the  probability  of  a  given  track,  say  the  first  track,  being  correctly  associated  with  a  mea¬ 
surement  when  conditioned  by  the  position  of  the  track,  i.e.,  Prob.(2?12|j)i).  For  the  sake 
of  simplicity,  we  assume  that  the  spatial  variation  matrix  V  is  diagonal  as  V  =  E 2I  and 
that  the  average  prediction  error  variance  matrix  S  is  (at  least  approximately)  diagonal  as 
5  =  a2I .  With  these  simplifying  assumptions,  the  conditional  probability,  Prob.(£12]yi),  is 
only  a  function  of  |jyj||.  Thus  the  purpose  of  this  appendix  is  to  calculate  the  conditional 
probability,  Prob.(£i2|||yi  ||  =  r). 

It  follows  from  (2.7)  and  the  above  assumptions  that 


Prob.(E„  |  Il9.ll  =r)  =  f  (erf  Aft II  =  r)  =  (-^)  f(p)dp  (C.l) 

where  /  is  the  probability  density  function  of  jjy:  —  y2|j  given  the  condition  ||yi||  =  r,  i.e., 

m  =  £«p  ’)  '•  (£)  <c-2) 

with  Iq(-)  being  the  modified  Bessel’s  function  defined  by 


00  i2" 


/0(l) 


It  follows  from  (C.l)  that 


P{Eu  I  ||yt||  =  r)  = 


r-  bt?)  £e*p  ( (^r)) u  (§)  ^ 


=  exp 


r2  \  r°°  2o2 

cV/o  E^ 


2E2 

(~r'\ 


erf(z)zexp  (  ~x2  }  I0 


ar 


E2 


-x  )  dx 


=  exp 


exp  (2^)  Qn?oCn^°°erf(a:)i2n+lexp  Hai0 dx 
00 


2\  00 


2E2 


(C.3) 


(C.4) 


n=0 
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2a-2  7 

where  a  =  — ,  Cn  = 


2n 


l  = 


s/l  dr 


,  and 


22n(n!)2 1 

Kn  d=  erf(x)z2n+1exp  ^az2j  dx 

The  sequence  {Kn}  can  be  calculated  recursively  as  shown  below:  First  we  should  note 

(2n-l)"J - 


f  exp  (—ax2)x2ndx  = 
Jo 


2  n+1 


(C.5) 


(C.6) 


Integrating  (C.5)  by  parts  and  using  the  above  equality,  we  have 


„2  n 


Kn  =  ^—Kn.x  + 

2na 


(2n  -  1)!! 


(a  +  l)n+1/2  22n+1(n!)2 


(C.7) 


where  (2n  —  1)!!  (2n  —  1)  •  (2n  —  3)  •  -3  •  1,  and 


K0  =  l  1  +  -=  (C.8) 

2  VS2  +  2<j2J  V 

We  can  calculate  the  conditional  probability  (C.4),  first  by  calculating  the  sequence 
{.fiTn}  by  (C.7)  and  then  by  putting  the  results  into  (C.4).  However,  an  approximate  but 
more  concise  expression  is  desirable.  In  order  to  obtain  such  an  approximate  expression,  we 
will  look  for  appropriate  upper  and  lower  bounds  on  the  series  Kn.  First  let  us  rewrite 
(C.7)  as 

a#2  TVn  -L 

2> 


K=jLK  i  72n  f(n  +  |) 

"  2na  2n+1(a  +  l)»+i/*  [r(n  +  l)]2 


Then  define  two  sequences,  {K'n}  and  { K"n },  by 


and 


where  p  =  — q  = 


K’  =  +  L 

n  n! 


K"  =  V-K"n- , 

n 


(C.9) 

(C.10) 

(C.U) 


2a  ’  ^  2(1  + a) 


,  and  K'o  =  K"o  =  K0.  Obviously,  we  have, 


£*“.<£  X.  <£  *'„ 


(0.12) 


n=0 


n^O 


n=0 


On  the  other  hand,  we  have  )T^L0  K"n  =  Koep,  and 


X>'„  =  /i:oep+l(ep-l)  +  ^(ep-l-p)  +  ... 

—  P  P 


n=0 


<  K0ep  +  (ep  -  1)  (-  +  ^  +  .. 

\P  P2  P3  > 
=  K0ep  +  a(ep  -  1) 


(C.13) 
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Therefore,  from  (C.12)  and  (C.13),  we  have. 


K0{t?  -l)<'£tKn<(K0  +  a)(e?  -  1) 


(C.14) 


Then,  assuming  j  «  S,  or  equivalently  a  <C  1,  we  can  make  the  approximations,  a  «  0, 
K0  ^  1,  and  hence 


£  * "  ~  ^  - 1 
n=l 

It  follows  from  (C.4)  and  from  this  approximation  that 


p{En  i  ||j/i|]  =  r)  a  exp  (  —  I  [K0  +  exp 


(AH—id) 


(C.15) 


£)’expfe) (C16) 


When  the  total  number,  N,  of  targets  is  random  and  has  a  Poisson  distribution  with 
mean  v ,  the  overall  correct  association  probability  can  be  calculated  as 


Prob.(jEe  |  ||jh||  =  r)  =  £Prob.(£e  |  ||^||  =  r>N  =  n)Prob.{JV  =  n} 

=  exp  (-1^(1  -  Prob.(£12  |  ||yx||  =  r))) 

=  expH(l)’'xp(i?)) 


(C.17) 


where  (3(r )  is  the  object  density  at  points  with  distance  r  from  the  object  group’s  center  and 
expressed  as 


/?(r)=2^CXp  [2& 


(C.18) 
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D.  Effects  of  False  Alarms 

In  this  appendix,  we  will  calculate  the  probability  of  correct  association  in  the  presence 
of  false  alarms.  As  we  did  in  Section  2.2  and  Appendix  B,  let  us  consider  an  object,  say 
the  first  object,  situated  at  origin,  in  the  sense  that  yx  =  0.  Then,  suppose  that  there  are 
N  —  1  other  objects  and  Nf  false  alarms,  both  N  and  Nf  being  random.  We  assume  that 
the  two  random  integers,  N  and  Nf ,  are  independent  and  have  Poisson  distributions  with 
mean,  u  and  v} ,  respectively.  Given  N  and  Nf,  the  joint  system,  {{y;}fl2,  {2/fA}i=i}  of 
other  objects,  {y,}^,  and  false  alarms,  {yf^4}^,  is  assumed  to  be  independent  and  has  an 
identical  uniform  distribution  on  the  disk  centered  >at  the  origin  with  radius  r. 

As  before,  let  Ec  be  the  event  in  which  the  first  object  is  associated  with  the  correct 
measurement.  Then,  since  the  objects  and  false  alarms  are  independent  from  each  other,  we 
have 

Prob.(£c)  =  Prob.(£rt)Prob.(£c/)  (D.l) 

where  Ect  is  the  event  in  which  there  is  no  transposition  among  objects  involving  the  first 
one,  and  Ec/  is  the  event  in  which  there  is  no  switch  between  the  true  measurement  from 
the  first  object  and  any  of  the  false  alarms.  Since  both  events  are  independent  from  each 
other,  the  results  of  Appendix  B  implies  ' 

Prob.(-Ert)  «  exp(— 7 r(3a2)  (D.2) 

1/  v  —  1 

where  (3  =  — -  « - —  is  the  object  density.  The  rest  of  this  appendix  describes  the  second 

7Ti  7TT* 

factor,  Prob^^y),  in  (D.l). 

Let  Pf  i  be  the  probability  that,  given  a  false  alarm  and  the  correct  measurement,  the 
false  alarm  is  not  associated  with  the  first  track,  i.e.,  the  likelihood  of  associating  the  true 
measurement  is  higher  than  the  likelihood  of  associating  the  false  alarm.  Then,  under  the 
condition  where  there  are  n  false  alarms  and  just  one  track,  the  probability  that  the  track 
is  associated  with  the  correct  measurement  (i.e.,  with  none  of  the  false  alarms)  is  (P/i)n- 
Therefore,  when  we  take  expectation  over  the  distribution  of  the  number  of  false  alarms,  we 
have 

ProME.,)  =  WO'  =  exp(-,,(l  -  P„))  (D.3) 

n=0  U- 

We  can  write  the  probability  Pf\  explicitly  as 

Pn  =  Prob.{||y!  -  yi||*_,  <  \\yf  -  y^.,}  (D.4) 

where  yx  is  the  true  measurement  which  should  be  associated  with  the  track  which  has  the 
prediction  at  yly  is  a  false  alarm,  and  5  is  the  average  prediction  error  variance  matrix 
satisfying  (B.4).  Moreover,  for  the  sake  of  simplicity,  we  assume  that  the  average  prediction 
error  variance  matrix  5  is  diagonal  as  S  =  a2I  at  least  approximately,  with  cr  being  the 
average  prediction  error  standard  deviation.  With  the  assumption  yx  =  0,  we  have 

Pn  =  Prob.{||yi ||2  <  M2}  =  ?iob.{\\yi\\  <  ||y/||}  (D.5) 
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Let  pi  =  j|yi||  and  pf  =  \\yj\\.  It  can  be  easily  shown  that  the  probability  density  functions 
of  pi  and  pf  are 

Pi(pi)  =  (- (D.6) 

and 

Tt(fl)  =  2jf  (D.7) 

r 

respectively,  with  the  same  support,  [0,r].  Therefore,  it  follows  from  (D.5)  that 


Pi{pi)P/{pf)dpidpf 


°n  =  J 

=  Jo  (J0'Pl(Tl)drx)pf(rridrf 


2<r*  ( 

=  1  “  7T  U-exp 


(-a) 


1  - 


2a 


(D.8) 


In  the  last  step,  an  approximation  based  on  R  a  has  been  used. 
Now,  it  follows  from  (D.3)  and  (D.8)  that 


Prob.(25cy)  =  exp 


=  exp  [~2Ff3FAa2) 


(D.9) 


where  fiFA  =  — r  is  the  false  alarm  density.  Therefore,  from  (D.l),  (D.2),  and  (D.9),  the 

7rr2 

overall  correct  association  probability  becomes, 

Prob.(£e)  =  exp  (-*(/?  +  20fa)&2)  (D.10) 
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E.  Orbital  Dynamics  and  Observation 


This  appendix  describes  the  dynamics  and  the  observation  functions,  /  and  h  in  Section 
3,  for  different  coordinate  systems  shown  in  Table  3-1. 

In  the  cartesian  coordinate  system,  eqn.  (3.1)  is  Newton’s  equation  for  the  two-body 
problem  with  x(t)  being  the  3D  position  and  velocity  at  time  t.  Thus  we  have 


/(* ) 


r 

\ 

f  »  1 

U 

t 

= 

V 

) 

- 1 

tn 

IT 

1 

_ 1 

(E.l) 


where  7  is  the  earth  gravity  constant.  The  observation  function  h(-,t)  at  time  t  maps  the 
vector  x(t)  —  x,(t)  into  the  angular  projection,  i.e.,  the  azimuth-elevation  pair,  where  x,(t) 
is  the  3D  position  of  the  sensor  at  the  time.  The  partial  derivatives  are  then  calculated  as 


df,  ,  df 

=  Tz 


u 

V 


uu‘ 


13^  3|IU|I2^ 


I 

0 


and1 


dh 

dx 


(x,t)=  ||z(i)  -  r.(f)|| 


-1 


e1  0 
QT  0 


(E.2) 


(E.3) 


where  {e,  3}  is  the  system  of  ortho-normal  vectors  which  are  orthogonal  to  x(t)  —  x,(t). 

When  we  choose  a  coordinate  system  such  that  z(t)  is  the  range  R(t ),  the  azimuth  A(t) 
and  the  elevation  E(t),  and  their  rates,  R(t ),  A(t )  and  E(t),  the  dynamics  represented  by  / 
becomes  very  complicated.  According  to  [16],  we  have 


R  =  R{E2  +  A2  cos 2  E)  -  - 
R  • 

A  =  — 2— A  -f-  2 AE  tan  E 
R 


7  sin  E 


R \ 


R2,{R  +  RsinE) 


R ^  sin  E 


-  1 


(E.4) 


-  R  • 

A2  .  r. 

7  cos  E 

'{R\3 

E  =  -2-E- 

—  sin  2  E  — 
2 

RR2, 

\rt)  1 

where  Rt  =  (R2  +  R\  +  2 RR,  sin  E)1?2  is  the  distance  of  the  object  from  the  earth’s  center, 
R,  is  the  distance  between  the  sensor  and  the  earth’s  center,  and  Rt  is  the  radius  of  the  earth. 
The  partial  derivatives  of  /  is  thus  very  complicated.  On  the  other  hand,  the  observation 
function  h  becomes  a  simple  projection  with  partials 


dh,  , 


010000 

001000 


(E.5) 


When  the  six  orbital  elements  are  used  for  the  object  state  variable  i(f),  the  dynamics 
become  very  simple.  For  example,  if  we  use  the  mean  anomaly,  the  mean  anomaly  rate,  the 


lI  is  the  identity  matrix  with  appropriate  dimension,  in  this  case  three. 
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angular  momentum  (or  the  eccentricity),  and  three  parameters  to  determine  the  orbital  plane 
and  the  perigee  direction,  ail  the  variables  except  one  (the  mean  anomaly)  are  constant,  and 
the  dynamics  are  simply  expressed  as  that  the  time  derivative  of  the  mean  anomaly  is  the 
constant  mean  anomaly  rate.  Thus  the  partials  of  /  are 

0  1  0  0  0  0 

0  0  0  0  0  0 

0  0  0  0  0  0 

0  0  0  0  0  0 

0  0  0  0  0  0 

0  0  0  0  0  0 

The  observation  function  h  becomes  the  composition  of  the  transformation  from  the  or¬ 
bital  elements  to  the  cartesian  and  then  from  cartesian  to  the  angular  measurement.  The 
partial  derivatives  can  be  obtained  by  post-multiplying  the  partials  of  the  orbital-element- 
to-cartesian  transformation  from  the  right  of  (E.3). 
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F.  Calculation  of  Probability  of  Measurement  Merging 


This  appendix  describes  the  derivation  of  (4.8),  i.e.,  the  probability  of  an  object  being 
unresolved  with  another  object  on  a  given  focal  plane  of  a  sensor,  in  Section  4.2.  For  the 
derivation,  we  will  use  the  following  mathematical  model: 


1. 


There  are  N  objects  on  the  focal  plane  where  N  is  a  random  integer  having  a  Poisson 
distribution  with  mean  v. 


2.  Given  the  number  N  of  objects,  the  objects  are  model  by  an  i.i.d.  system, 

with  a  common  distribution  which  is  uniform  on  the  disk  centered  at  the  origin  and 
with  radius  r. 


3.  The  probability  of  merging  two  objects,  at  yx  and  y2,  are  given  by  pm(||yi  —  y2||),  where 
pm  is  a  continuously  decreasing  function  such  that  pm(0)  =  1  and  pm{p)  =  0  for  p  >  8 
with  S  being  the  sensor  resolution.  In  particular,  we  will  adopt  a  piecewise  linear  model 
described  by  (4.6)  in  Section  4.2.  (See  Figure  4-1.) 

4.  We  assume  that  v  1  and  r  >•  6  but  the  object  density 


is  comparable  to  1  with  respect  to  the  resolution  cell  area,  62,  i.e., 


is  comparable  to  1. 


Let  Em  be  the  event  in  which  an  object,  say  the  first  object,  zx,  is  merged  with  (unre¬ 
solved  from)  at  least  one  of  other  objects,  z2,...,zjv,  i.e., 


Em  = 


(F.l) 


where  E'm  is  the  event  in  which  the  first  object  is  merged  with  the  i-th  object.  Unfortu¬ 
nately,  it  is  not  easy  to  calculate  the  probability,  Pm,  of  the  event  Em  mainly  because  the 
Em'^’s  are  not  disjoint.  In  order  to  overcome  this  difficulty,  we  will  approximate  the  event 
Em  by  Em ■  where  i  is  the  index  of  the  nearest  neighbor  object  zj  from  zx.  The  validity  of 
this  approximation  is  tested  by  Monte  Carlo  runs,  the  results  of  which  are  shown  in  Section 
4.2.  (See  Figure  4-2.)  With  this  approximation,  we  may  calculate  the  probability,  Pm,  of 
measurement  merging  as 


Pm 


def 


Prob.(Fjvr) 

Prob.f  EfrTl) 

Prob .{E^  |  N  >  2)Prob .{N  >  2} 

I  Pm{p)q{p)dp 
Jos 

I  Pm(p)q(p)dp 

Jo 


(F.2) 
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i 


where  q  is  the  probability  density  function  for  the  distance  to  the  nearest  neighbor  defined 
by 

q{p)dp  Prob.{||2i|j  6  dp  |  N  >  2}  (F.3) 

In  the  last  step  of  (F.2),  we  used  the  fact  u  1. 

Thus  we  will  first  calculate  the  distribution  of  the  distance  from  the  first  object  Z\  to 
its  nearest  neighbor  rj.  Assuming  a  large  enough  radius  r  of  the  extent  of  the  objects,  we 
can  let  z\  =  0.  Other  objects,  Z2,...,z.n,  form  an  i.i.d.  system  with  the  uniform  distribution 
on  a  disk,  and  hence, 

J  r  n 

(F.4) 


p(p)dp  d=  Prob.{|jz,||  €  dp}  =  2^-dp 


for  p  6  [0,r].  Therefore,  the  distribution  of  the  distance  to  the  nearest  neighbor  object  (from 
the  first  object)  can  be  expressed  as 


pn(p)dp  ^  Prob.{||z?|!  €  dp  |  N  =  n  +  1}  =  2n?¥ -/T~)  dp 


(F.5) 


for  p  6  [O.r]  and  for  n  >  1.  Since  we  are  assuming  that  the  number  N  of  object  has  the 
Poisson  distribution  with  mean  v,  we  have 


i/n+1)  2 np(r2  -  p2)n~'1 


<l(p)  = 


ll 

3 

+  1)1  r2n 

“  ^(n+l) 

5(n+  !)! 

_  2  P  - 

K-s)-: 

e(1-£)l'  + 1 

r2 

[ev  -  1  -  v)  | 

:-s>‘ 

(F-6) 


C  P  exp  (-tt/V) 

for  p  €  [0,5],  where  C  is  the  appropriate  normalizing  constant.  In  the  last  step  in  the  above 
equation,  we  have  used  approximations  based  on  assumptions  v  1  and  r  6.  Figure  F-l 
shows  the  probability  density  function  q(-)  for  several  different  values  for  /?. 

It  follows  from  the  two-object  merging  model,  pm,  defined  by  (4.6),  and  (F.2)  to  (F.6) 

that 

Pm  «  /  Pm{p)q(p)dp 

Jo 


rS/2  ,  r 

2 *0  J  pe~*0>  dp  +  2J* 


*  6~P---*0S 


i/2  6 


-pe 


dp 


(F.7) 


=  1  - 


erf  (^2 -  erf 


which  proves  (4.8). 

We  should  note  here  that  the  probability  of  measurement  merging  (or  unresolved  mea¬ 
surements),  Pm,  is  a  function  of  the  object  density  f3,  or  the  object  density  normalized  by 
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the  sensor  resolution,  i.e.,  (382.  As  expected,  we  have  Pm  T  1  35  T  °°  and  Pm  1  0  as 
( 36 2  1  0.  We  have  assumed  the  object  density  j3  to  be  uniform.  If  the  density  is  not  uniform 
(as  discussed  in  Appendix  C),  the  calculation  of  the  probability  density  functions,  (F.5)  and 
(F.6),  may  become  very  complicated,  and  hence,  we  may  not  have  a  concise  expression  such 
as  (4.8).  However,  if  the  change  of  the  object  density  is  gradual  enough  relative  to  the  sensor 
resolution  8 ,  the  probability  Pm  may  be  accurate  when  0  is  replaced  by  the  local  object 
density,  and  the  overall  (average)  probability  may  be  given  by  appropriate  integration. 
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